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Abstract 



In this work we investigate the quantum dynamics of three different configu- 
rations of nonhnear cavity systems. We begin by carrying out a quantum analysis 
of a dc superconducting quantum interference device (SQUID) mechanical displace- 
ment detector comprising a SQUID with a mechanically compliant loop segment. 
The SQUID is approximated by a nonlinear current-dependent inductor, inducing an 
external flux tunable nonlinear Duffing term in the cavity equation of motion. Ex- 
pressions are derived for the detector signal and noise response where it is found that 
a soft-spring Duffing self-interaction enables a closer approach to the displacement 
detection standard quantum limit, as well as cooling closer to the ground state. Next, 
we consider the use of a superconducting transmission line formed from an array of 
dc-SQUIDs for investigating analogue Hawking radiation. We will show that bias- 
ing the array with a space-time varying flux modifies the propagation velocity of the 
transmission line, leading to an effective metric with a horizon. As a fundamentally 
quantum mechanical device, this setup allows for investigations of quantum effects 
such as backreaction and analogue space-time fluctuations on the Hawking process. 
Finally, we investigate a quantum parametric amplifier with dynamical pump mode, 
viewed as a zero-dimensional model of Hawking radiation from an evaporating black 
hole. The conditions are derived under which the spectrum of particles generated 
from vacuum fluctuations deviates from the thermal spectrum predicted for the con- 

ii 



ventional parametric amplifier. We find that significant deviation occurs once the 
pump mode (black hole) has released nearly half of its initial energy in the signal 
(Hawking radiation) and idler (in-falling particle) modes. As a model of black hole 
dynamics, this finding lends support to the view that late-time Hawking radiation 
contains information about the quantum state of the black hole and is entangled with 
the black hole's quantum gravitational degrees of freedom. 
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Introduction and organization 

1.1 Organization of this thesis 

This thesis is naturally organized into three chapters that are relatively indepen- 
dent from each other: 

• Displacement detection and cooling of a mechanical oscillator using a nonlinear 
microwave cavity detector (Chapter |2| 

• Analogue Hawking radiation in a dc-SQUID array transmission line (Chapter|3]) 

• The Trilinear Hamiltonian: Modeling Hawking radiation from a quantized source 
(Chapter g 

Given the self-contained nature of the topics, each chapter includes its own in- 
troduction and conclusion sections. However, the work in the later chapters relies 
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upon the knowledge gained in the previous chapters and is the underlying reason for 
the ordering of research topics presented here. To highlight these relations, Chapters 
[3] and |4] contain "Motivation" sections which attempt to elucidate the connections 
between the various research topics. Given the overall theme of nonlinear cavity 
systems, in the next section we briefly introduce the field of superconducting circuit 
cavity systems that provides the foundation for the work presented in this thesis. 

1.2 Prologue 

Since the early days of quantum mechanics the interaction between two-level sys- 
tems and harmonic oscillators has played a central role in our understanding of physics 
below the macroscopic level. Beginning with the thought experiments of Einstein and 
Bohr, scientists have devoted significant effort towards understanding the dynamics 
of these most basic elements of quantum theory. Unfortunately, as is often the case 
with physics, it would be another fifty years before these conceptual ideas could be 
brought to reality. It wasn't until the 1970's and the advent of Cavity Quantum 
Electrodynamics (CQED) that physicists could fabricate and control the interaction 
between an atom (two-level system) and an optical oscillator (harmonic oscillator) 
while at the same time removing unwanted environmental effects. Since then, CQED 
has contributed greatly to our fundamental understanding of the interaction of mat- 
ter with quantized electromagnetic fields, the physics open quantum systems, and 
decoherence (see Ref. [1] and references therein). Additionally, these same systems 
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will undoubtably be utilized in the future realization of a quantum computer. The 
experiments made possible by CQED have brought our understanding full circle to 
where it is possible to test the very foundations upon which quantum mechanics rests. 

Building upon the success of CQED, this last decade has seen a large effort de- 
voted to exhibiting complementary effects in solid-state superconducting circuits, a 
field of research known as circuit-QED[21 E]- While similar to CQED, the use of 
superconducting circuits to fabricate artificial atoms and cavities has several unique 
advantages. To begin, the effective two-dimensional coplanar microwave transmission- 
line cavity^ geometries used in superconducting circuits have mode volumes that are 
markedly smaller than the corresponding three-dimensional Fabry-Perot cavities used 
in their optical counterparts. This modal volume is significantly less than the cubic 
wavelength corresponding to a microwave photon in vacuum (A = 30 cm at 1 GHz) 
and consequently leads to a greatly enhanced electric field strength [3J and corre- 
spondingly large amplitude vacuum fluctuations [SJ E]- Additionally, the small size 
and low-temperature environment of these resonators has routinely produced quality 
factors of Q ~ 10^ or higher |1] making it possible to reach the ultra-strong coupling 
regime [3 |B]. 

The artificial atoms (qubits) to which the cavity resonator interacts are composed 
of highly nonlinear Josephson junction circuit elements[9] that sufficiently modify the 
harmonic potential of a typical LC-circuit allowing for the isolation of the lower two 
energy levels p^. The energy spacing between levels can be fabricated to a desired 



3 



Chapter 1: Introduction and organization 



value and then modified in situ by the apphcation of an external magnetic flux bias. 
At present there are a variety of experimental realizations of superconducting atoms 
based on charge, phase, and flux qubits, each distinguished by the degree of freedom 
used in generating the approximate two-level potential|llt Since the original 

demonstration of a solid-state qubit [,13J the coherence times of these devices have 
gone from a few tenths of nanosecond up to ~ one microsecond [12], representing 
an increase in four orders of magnitude over the last ten years. This rapid rise 
in coherence times has helped to positioned solid-state qubits and the circuit-QED 
architecture at the forefront of the race for a scalable quantum computer. 

To date, circuit-QED has resulted in unparalleled access to the coherent interac- 
tions between qubit and cavity resonator. It is now possible to control and observe 
photon number states in a microwave resonator using a coupled qubit [T ^ [T5 | [TB]. gen- 
erate and perform tomography of arbitrary superposition states[l7], and measure the 
decay of these states due to the ever present environment [IH] . Qubits coupled to cav- 
ity resonators [19] can find applications in single- atom lasing[20] and fluorescence [21] , 
quantum amplifiers [22] . coupled qubit dynamics, and the demonstration of quantum 
computing gates and algorithms (an exhaustive list of references on qubit coupling 
and computation can be found in Ref. [T2]). 

In this thesis we will make use of the circuit-QED architecture in exploring the 
detection and cooling of a mechanical oscillator by a cavity detector, and as analogue 
models for investigating classical and quantum gravitational physics. The prepara- 
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tion and manipulation of quantum states of a mechanical resonator has many po- 
tential applications in fundamental physics and applied sciences [231 [2H 125] . This 
has stimulated many groups to focus on coupling mechanical resonators to both 
superconducting[2Sl [23 ESI ESI EDI ED E21 E3] and optical cavities|Ml ES] with the 
hope of reaching the quantum ground state in the mechanical element, a goal that 
has recently been achieved However, the push to see quantum effects in pro- 
gressively larger systems requires ever more advanced methods to reach this goal. In 
Chapter |2] we will demonstrate the detection and cooling of a mechanical resonator 
coupled to a nonlinear microwave cavity. In Chapters [3] and |4] we will focus on simu- 
lating the gravitational physics of black holes using superconducting circuits. The use 
of superconducting elements in the investigation of analogue gravitational effects is a 
previously unexplored area that allows the possibility of observing otherwise elusive 
quantum gravitational phenomena in the laboratory. 
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Displacement detection and 
cooling of a nanomechanical 
resonator using a nonlinear 
microwave cavity detector 

2.1 Introduction 

Recently there has been interest in exploiting the nonlinear dynamics of nano- 
electromechanical systems (NEMS) for amplification. |37l l38t [39] The use of nonlinear 
mechanical resonators to some extent parallels investigations with systems comprising 
purely electronic degrees of freedom, such as nonlinear superconducting devices in- 
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corporating Josephson Junctions (J J). [101 El 112] For example, it was shown that the 
bistable response of an RF-driven JJ can be employed as a low noise, high-sensitivity 
amplifier for superconducting qubits. [40j A similar setup consisting of a J J embedded 
in a microwave cavity was used to measure the states of a quantronium qubit.[15] 
where the relevant cavity mode was found to obey the Duffing oscillator equation. 

One area of nanomechanics that has yet to fully explore the possibility of ex- 
ploiting nonlinearities for sensitive detection involves setups in which a nanome- 
chanical resonator couples either capacitively[26l |29l [30] or inductively[27l [2Ej to 
a superconducting microwave transmission line resonator, combining elements from 
both the above-described NEMS and superconducting systems. Such setups are in 
some sense the solid-state analogues of optomechanical systems, which ponderomo- 
tively couple a movable mirror to the optical field inside a cavity using radiation 
pressure. [l5l HHl HZl HHl HHl [50l EH |52l [53] In both areas, the focus has primarily been 
on operating in the regime where the cavity and resonator behave to a good approxi- 
mation as harmonic oscillators interacting via their mutual ponderomotive coupling. 
However, in the case of microwave cavities, introducing an embedded JJ,[13] or sim- 
ply driving the cavity close to the superconducting critical temperature. [5^ results in 
the breakdown of the harmonic mode approximation; nonlinear dynamical behavior 
of the cavity must be taken into account. Furthermore, the ponderomotive coupling 
term between the microwave or optical cavity mode and mechanical mode is by itself 
capable of inducing strong, effective nonlinearities in the respective mode equations. 
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In optical systems, such nonlinearities can manifest themselves in the appearance of 
a bistable (or even multistable) region for the movable mirror. By restricting 

ourselves to linear microwave cavities, we are overlooking a range of nonlinear phe- 
nomena that might enable a closer approach to quantum-limited detection, as well 
as cooling of the mechanical oscillator closer to its ground state. As an illustration, 
consider the phase sensitive Josephson parametric amplifier. [57t [58 | l59 | 160] which ex- 
ploits the nonlinear effective inductance of the JJ to perform (in principle) noiseless 
amplification and quantum squeezing of the respective complimentary quadrature 
amplitudes of the signal oscillator. 

In this chapter, we will go beyond the usually considered ponderomotively-coupled 
two oscillator system to include a Duffing nonlinearity in the microwave cavity mode 
equations. The closed system model Hamiltonian describing the nonlinear microwave- 



coupled mechanical oscillators is given by Eq. (2.26). The nonlinear microwave mode 
is externally driven with a pump frequency Up that can be detuned from the transmis- 
sion line mode frequency ujt- Our investigation will focus on the nonlinear amplifier 
created by embedding a dc-SQUID displacement detector into a superconducting mi- 
crowave transmission line. [27] This has the advantage of significantly larger coupling 
strengths[n] as compared with existing geometrical coupling schemes. [2H1 ESI ED] The 
displacement detector comprises a SQUID with one segment consisting of a doubly- 



clamped mechanical resonator as shown in Fig. 2.1 The net fiux, and therefore 



circulating current, is modulated by the mechanical motion, providing displacement 
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transduction. The capacitively-coupled pump/probe feedline both drives and pro- 
vides readout of the relevant transmission hne resonator mode amphtude (or phase). 
We will assume transmission line losses are predominantly due to coupling with the 
feedline, and that the pump drive is coherent. The main irreducible noise source 
is therefore microwave photon shot noise from the drive that acts back on the me- 
chanical oscillator via the intermediate nonlinear microwave resonator and SQUID. 
Environmental influences on the mechanical oscillator other than that due to the 
SQUID detector are simply modelled as a free oscillator thermal bath. By operating 
the amplifier well below the superconducting critical temperature, and with trans- 
mission line currents less than the SQUID JJ's critical current threshold, resistive 
tunneling of electrons and the associated noise is a negligible contribution. Similar 
setups involving J J elements have been considered previously. [271 [28t |6T| [62] 

With JJ plasma frequencies assumed to be larger than both the mechanical and 
transmission line fundamental mode frequencies, the SQUID can be considered as a 
passive, effective inductance element that depends on both the externally apphed flux 
and drive current. The effective inductance can therefore be freely tuned by varying 
these external parameters. Previously, we considered only the lowest, zeroth order 
expansion of the inductance with respect to the current entering (or exiting) the 
SQUID, [27] yielding the usual ponderomotively-coupled double harmonic oscillator 
system. In this companion investigation, we include the next leading, quadratic 
order term, resulting in a nonlinear current dependent inductance. Provided that 
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the current magnitude is small as compared with the JJ's critical current, neglecting 
higher order terms should not introduce significant errors. The nonlinear inductance 
induces an effective Duffing (i.e., cubic) self-interaction term in the microwave mode 
equations of motion. The results presented here apply to a broad class of bosonic 
detector, which includes optomechanical amplifiers with nonlinear cavities [63j that 



are describable by Hamiltonian (2.26). A related analysis of quantum noise in a 



Duffing oscillator amplifier is given in Ref. 



The chapter is organized as follows. In Sec. 2.2 we first derive the truncated 



Hamiltonian (2.26) that describes the closed system dynamics of the coupled cavity 
and mechanical resonator fundamental modes. We then derive the quantum Langevin 
equations of motion that describe the open system dynamics in the presence of the 



we 



pump/probe line and mechanical oscillator's external environment. In Sec. 2.3 
find expressions for the detector signal response and noise using a semiclassical treat- 
ment of the detector's linear response to the external noise input signal driving the 
mechanical oscillator. Section 12.41 determines the critical drive current for the on- 



set of bistability (not to be confused with the J J critical current). Sections 2.5 and 



2.6| discuss the effects of the microwave mode Duffing nonlinearity on mechanical 
mode displacement detection and cooling, respectively, giving illustrative examples 
assuming achievable device parameters. Section [2^1 briefiy concludes, while the more 



technical aspects of the signal and noise term derivations are relegated to Appendix [Xj 



Source code for the numerical analysis in Sees. 2.5 and 2.6 is given in Appendix B.l 
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2.2 Equations of Motion 



2.2.1 Closed System Hamiltonian 



The displacement detector scheme is shown in Fig. 2.1 The device consists 
of a striphne resonator (transmission hne) of length / bisected by a SQUID. The 
transmission line is characterized by an inductance and capacitance per unit length 
Lt and Ct respectively. The SQUID comprises two identical Josephson junctions with 
critical current Ic and capacitance Cj. A segment of the SQUID loop is mechanically 
compliant, forming a doubly clamped resonator of length /qsc- A similar setup without 
the microwave cavity has recently been constructed [65]. We only take into account 
mechanical fundamental mode displacements in the plane of the loop and assume 
that the resonator can be modeled effectively as a harmonic oscillator with the y 
coordinate giving the center of mass displacement. The magnetic flux threading the 
loop is given by $cxt(y) = ^c^t{0)+\Bc^Joscy, where $ext(0) = $ext is the flux with the 
mechanical oscillator fixed at y = 0, i?ext is the externally applied field in the vicinity 
of the resonator, and A is a geometrical factor that corrects for the non-uniform 
displacement of the oscillator along its length. The coupling between mechanical 
oscillator and external heat bath is characterized by the oscillator amplitude damping 
rate jbm, while the pump-probe line-transmission line coupling is characterized by the 
transmission line amplitude damping rate jpT- In what follows, we will assume weak 
couplings (i.e., large quality factors for the transmission line and mechanical oscillator) 
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and also that the dominant dissipation mechanism for the transmission hne is due to 
its couphng to the pump-probe hne, 7^^- 




-l/2< *-l/2 

Figure 2.1: Layout of the dc SQUID displacement detector. The dimensions of the 
SQUID have been enlarged relative to the transmission line to show the key charac- 
teristics employed in the analysis. 



In analyzing the SQUID dynamics, an appropriate choice of variables is 7-1- = 
(01 ± 02)/2, where 0i and 02 are the gauge invariant phases across the Josephson 
junctions. [66] while for the transmission line we choose the phase field 0(x, t). The 
transmission line current and voltage are described in terms of 0(a;, t) using the stan- 
dard telegraphic relations: 

Vi., ^ ,2.) 

where $0 = h/{2e) is the flux quantum. Assuming that the SQUID can be lumped 
at the midpoint a; = of the transmission line, the equations of motion for the closed 
system comprising the SQUID, transmission line and mechanical oscillator are given 
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by[67] 



92 



1 92 



U j 7_ + cos (7+) sin(7_) + 2/3 



L 



7_ — TT n + 



$0 



and 



7+ + sm(7+) cos(7-) + 1^^^ = 0, 



my + mul^y - -%A5ext^osc7- = 0, 



(2.3) 
(2.4) 
(2.5) 

(2.6) 



where ojj = y27r Jc/ {C j$o) is the plasma frequency of the Josephson junctions, /3l = 
27rL/c/$o is a dimensionless parameter with L the SQUID loop self-inductance and 
Ic the Josephson junction critical current, and where n takes on integer values arising 
from the requirement that the phase around the loop be single- valued. Equation 



(2.3) is the wave equation for the transmission line, equation (2.4) describes the 



current circling the loop, which depends on the external flux and oscillator position. 



equation (2.5) describes the average current in the loop, and Eq. (2.6) is Newton's 



second law for the mechanical oscillator with Lorentz force acting on the oscillator. 
The current and voltage across the SQUID must also obey the boundary conditions 



90(±//2,t) 



0; 



90(0", t) 90(0+, t) 



7+ - 



dx dx dx 

L 9^0(0, t) 90(0", t) 90(0+, t) 



ALt dtdx 



dt 



dt 



(2.7) 
(2.8) 



Using Eqs. (2.3)-(2.8), we shall now derive approximate equations of motion de- 



scribing a single mode of the transmission line interacting with the mechanical os- 
cillator, where the form of the interaction between the two oscillators is governed 
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by the SQUID parameters and boundary conditions. Assume that the following 
conditions are satisfied: (a) ujj ^ ujt ^ <^m- (b) (3l ^ 1- (c) |/(0,t)//c| ^ 1. 
(d) lA-Bext^oscZz/'^'ol ^ 1- Condition (a) states that the SQUID plasma frequency is 
much larger than the transmission line mode frequency of interest, ut, and conse- 



quently we shall ignore the SQUID inertia terms in (2.4) and (2.5). Condition (b) 
allows us to neglect the SQUID loop self inductance and, together with (a), elimi- 
nate 7-1- from the equations by expressing them in terms of the transmission line and 
oscillator coordinates as series expansions in f3L- Conditions (c) and (d) allow us to 
expand the above equations in the transmission line current l{0,t) = I{t) at x = 
and in the oscillator displacement y. Keeping terms to first order in y and to leading. 



second order in /, Eq. (2.6) for the mechanical oscillator becomes approximately 



my + mul^y - LqJ'^ /2 = 0. 



(2.9) 



The voltage boundary condition (2.8) can be expressed as 



d(p{0-,t) 90(0+, t) 



dt 



dt 



(2.10) 



where L{I,y) is the effective inductance, which expanded to second order in / takes 
the form 

L{I,y) = Loo + L20 {I/Icf + Loiy, (2.11) 



14 



Chapter 2: Displacement detection and cooling of a nanomechanical resonator using 
a nonlinear microwave cavity detector 



where the Lij coefficients are defined as 







967r/c 



4/, 



sec (7r$ext/$o) 
sec^ (7r$ext/$o) 

sec (7r<l>ext/$o) tan (7r<l>ext/$o) 



(2.12) 
(2.13) 
(2.14) 



Note that we have neglected the Py term in (2.11), restricting ourselves to the lead- 



ing order coupling only between the transmission line and mechanical oscillator, as 



already stated. The above equations differ from those of the prequel [27j through the 
inclusion of the nonlinear, leading order current-dependent contribution [L2o{I/IcY] 
to the effective inductance L{I,y). 



The nonlinear voltage boundary condition (2.10 ) with inductance given by Eq. (2.11 ) 



generates frequency tripling harmonics of the transmission line resonator mode. Omit- 
ting for the time being the mechanical oscillator degree of freedom, a trial perturba- 



tive mode solution to the wave equation (2.3) that includes the leading harmonic and 



solves the current boundary conditions (2.7) is the following: 



(f){x,t) 



+A cos(wt + Lp) COS [k{x - 1/2)] + aA^ cos(3a;t + Sip) cos [3k{x - 1/2)] ; a; > 

-A cos(cut + ip) cos [k{x + 1/2)] - aA^ cos{3ujt + 3(p) cos [3k{x + 1/2)] ; a; < 0, 

(2.15) 



where k = k^^^ + k^-^^ and to = \k\/\fhxCT- The coefficients a, k'^-^^ and k^^^ are 



determined by substituting Eq. (2.15) into the voltage boundary condition (2.10) and 



solving perturbatively to order A? ^ with k^^^ scaling as . We obtain: a = —1/48, 



(A;(°)//2)tan (A;(°)//2) = 



(2.16) 
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and 

fc(i)/ = --C^A\k^%^sm^ (fc(°)//2) , (2.17) 



where 



Considering the transmission hne phase field at the location x = —1/2, where the field 



is pumped and probed (see Fig. 2.1), the perturbative solution (2.15) can be obtained 



from the following single mode equation for 0(— Z/2,t) = 0(t): 



"^'1 + + ^4 f 1 - 18C^ sin (fcW//2)] 03 



dt^ 12 

1 (l - 2C^ [(A:(°)//2) sin (A;(°)//2)] ^) ^ = 0, (2.19) 



where ut = |^*-°''|/V-^tC't- The awkward nonlinear term (p^ can be eliminated by 
redefining the phase mode coordinate as = '?/'(l + Tip"^), provided |r|0^ ^ 1, where 

r = -1 (l - 2C^ [{k^'H/2) sin (A:W//2)]') . (2.20) 



The mode equation (2.19) in terms of the redefined phase coordinate ip then becomes 



+ 4^ - ^ooie [{k^'^l/^) sin (fc(°)//2)] ' tf;^ = 0. (2.21) 

Thus, embedding a SQUID in a microwave transmission line induces a cubic nonlin- 
earity in the effective single mode equations (under the conditions of small currents 
as compared with the Josephson junction critical current), resulting in the familiar 
(undamped) Duffing oscillator. 
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We now restore the mechanical degree of freedom y{t) by assuming that for small 
and slow displacements [conditions (a) and (d) above], the interaction with ip can be 
obtained by expanding ujt [through its dependence on $ext(l/)] to first order in y in 



Eq. (2.21) to obtain 



i, + u:lij - -4^3^(A;W//2)sin(A;(°)//2)]'V'' 
3 



Ai? Josey ^0 (7,$,,,/$o) sec (7r$ext/$o) 4^. (2.22) 



Equation (2.9) for the mechanical oscillator, together with Eq. (2.22) for the phase 



coordinate, follow from the Lagrangian: 

^'<P,y,'^,y) = ^mf- ^mujl^y^ + hjrl f^j sin2(A;o//2) 



($o/7r) AttLtIIc 



uoli>'^{t) 

+ ^4C' [(A;^°)//2) sin (A;WV2)]' V-'} • (2.23) 

We now introduce the phase momentum coordinate p^ = dC/dip = m^ip and raising 
(lowering) operators 



^ - ^ (m^UTi^ =F iPi,) (2.24) 

= ^ (mUmyTipy) (2.25) 

satisfying the usual commutation relations, where the effective phase mass is = 
\CtI ($o/27r)' sin^ {kQl/2). In terms of the raising (lowering) operators, the Hamilto- 
nian operator is 

H = hhjTO^aT + ^hwrKdia:^ + arY + hwma:l;^am + ^HurKTmiaT + "'rYi^^m + «m), 

(2.26) 
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where, for notational convenience, hats on the operators and the minus superscript on 
the lowering operators will be suppressed from now on. The parameter characterizing 
the strength of the interaction between the transmission line mode and mechanical 
oscillator mode is 

KTm = ^%¥^r^ tan (vr$e.t/$o) sec (vr$e.t/$o) , (2.27) 

where A^p = ^Jfij {2mUm) is the zero-point displacement uncertainty. The parameter 
Kd characterizing the strength of the Duffing nonlinear term takes the form 

' '{2ef/{2CTl)' 



(2.28) 



which has been written in such a way as to make clear its various dependencies. 
In particular, Kd depends essentially on the cube of the ratio of the linear SQUID 
effective inductance Lqo to transmission line inductance L^Z, as well as on the ratio of 
the single Cooper pair charging energy to the microwave mode photon energy of the 
transmission line. Since the strength and sign of the linear SQUID inductance depends 



on the external flux $ext [see Eq. (2.12)], it is possible to vary the strength as well 
as the sign of the Duffing constant by tuning the external flux either side of $o/2. 
Thus, we can have either spring hardening or spring softening of the transmission 
line oscillator mode. Previously this flux tunability was observed in the readout of 
a persistent current qubit.[l2] Note, however, that the perturbative approximations 



that go into deriving the above Hamiltonian (2.26) do not allow too close an approach 



to the singular half- integer flux quantum point. In particular, the validity of the 
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expansions in Ic and /3l properly require the following conditions to hold: 

I 



sec (7r<l>ext/$o) 

J- r 



< 1 (2.29) 

|/3Lsec(7r$ext/$o)| < 1- (2.30) 



As already noted, Eq. (2.26) without the Duffing nonlinearity coincides with the 
Hamiltonian commonly used to describe the single mode of an optical cavity inter- 
acting with a mechanical mirror via the radiation pressure. However, we have just 
seen that embedding a SQUID within a microwave transmission line cavity induces 
a tunable Duffing self-interaction term as well; it is not so easy to achieve a similar, 
tunable nonlinearity in the optical cavity counterpart. 

2.2.2 Open System Dynamics 

Up until now we have considered the transmission line, SQUID and mechanical 
oscillator as an isolated system. It is straightforward to couple the transmission line 
to an external pump-probe feedline and mechanical oscillator to a thermal bath using 
the 'input-output' formalism of Gardiner and Collett.[68] Assuming weak system-bath 
couplings justify making the rotating wave approximation (RWA), and furthermore 
making a Markov approximation for the bath dynamics, the following Langevin equa- 
tions can be derived for the system mode operators in the Heisenberg picture: 



dt h\l 2mujm 



-iUJrnam + T\ TT" -^ext(^) " iUrKTmO'TO'T 



-76ma„^-zV2^e''^""4'^(t) (2.31^ 
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— ioJTCiT — ioOT^^dO'T^TO'T ~ i^TKTmO'T{cL^ + Cim) 



and 

dttT 
~dt 

-IpraT - i^/2^e'^''^af{t), (2.32) 

where jbrn is the mechanical oscillator amplitude damping rate due to coupling to 
the bath, 7^7- is the transmission line mode damping rate due to coupling to the 
pump-probe line, and we have also assumed that the small Duffing coupling and 
transmission line- mechanical oscillator coupling justify applying the RWA to the 
transmission line mode operator terms. The 'in' bath and probe line operators are 
defined as 

aj°(t) = / rfwe~^"(*-*°)ai(w,to), (2.33) 

\/2n J 



where t > to, with the states of the pump-probe line and oscillator bath assigned at 
to, interpreted as the initial time in the past before the measurement commences. For 
completeness, we have also included a classical, external time-dependent force -Fext(t) 
acting on the mechanical oscillator, although we shall not address the force detection 
sensitivity in the present work. 

It will be convenient to work with the Fourier transformed Langevin equations. 



With 0{tu) = J^^due'^^^Oit), Eqs. (2.31) and (2.32) become 



00 



+ ^^^ I dJ [aT{uj')a+{uj' + a+(cj')«T(c^' + w)] |^ (2.34) 
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and 

f^OO POO 

X dcu' duj"a:^{uj")aT{co')aTico + uj" - uj') 



J —oo J — oo 



27r 



doo'aTiuj') [omioo - oo') + a+{uj' -oo)]\ ■ (2.35) 



2.3 Detector Response 

The probe line observables are expressed in terms of the 'out' mode operator: 

a7*(t) = ^ / duje-'''^'-''^ap{uj,ti), (2.36) 



where ti > t. The 'out' and 'in' probe operators are related via the following useful 
identity: [68] 



a-^*(t) = -Zv^e-*<^--aT(t) + afit), (2.37) 

which allows us to obtain the expectation value of a given observable once ar(t) is 
determined. As an illustrative expectation value, we shall consider the variance in 
the probe line reflected current in a given bandwidth 6uj centered about the signal 
frequency of interest a;o:[27] 

1 /"^»+'5'^/2 dujidu2 / 2 sin [{ui - uj2)Tm/2] 



xi(a-*(a;0«r^(^2) + a-*+(a;2)a°^*(u;i)), (2.38) 



where, in addition to the ensemble average, there is also a time average denoted by 
the overbar, with the averaging time taken to be the duration of the measurement 
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Tm, assumed much longer than all other timescales associated with the detector dy- 
namics. In particular, time averaging is required when F^xtit) has a deterministic 
time dependence. [27] Expectation values of other observables, such as the reflected 



voltage variance and reflected power are simply obtained from Eq. (2.38) with appro- 



priate inclusions of the probe line impedance Zp = ^jLpjCp. P°"* = ([(5V°"*]2)/Zp 



From the form of the -ft'Tm coupling term in Eq. (2.32), we can see that the 
motion of the mechanical resonator modulates the transmission line frequency, and 
thus a complimentary way to transduce displacements besides measuring the cur- 
rent amplitude, is to measure the frequency- dependent, relative phase shift between 
the 'in' pump drive current and 'out' probe current using the homodyne detection 
procedure. [69] While we shall focus on amplitude detection, the homodyne method 
can be straighforwardly addressed and is expected to give similar results for the quan- 
tum limited detection sensitivity. 



Substituting Eq. (2.34) into (2.35), we obtained the following single equation for 



the transmission line mode operator a^: 



ar(w) = / dJ aTijjJ — iJ^Aiijj.iij') ^ \ du'B{uj,uj')aT{uj — u)') 

J —oo J — oo 

/oo 
dco" [aT(a;")4(w" - w') + a^{u")aT{uj" + to')] 
-oo 

poo poo 

+D{u) / du" / du'a^{u")aT{uj')aT{uj + tu" — uj') + C{u) ,{2.39) 
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where 



A{uj,uj') 



utK' 



Tm 



1 



27r u -ujt + ilpT 



1 1 



47r W - + ^7pT L^' - l^m + ilhn 

St{uj) 



-UJ' - Urn- ilbn 



C{UJ) 



(2.40) 

, (2.41) 
(2.42) 



and 



with mechanical signal operator 



27r tU — LOT + ilpT 



(2.43) 



:Fext(u;) 



(2.44) 



and noise operator 



St{u) = y2^e^^^-a;(a;) 



(2.45) 



For small signal strength, it is assumed that Eq. (2.39) can be solved as a series 



expansion up to first order in A{u, u'), giving the usual linear-response approximation. 
I.e., aT{u)) ~ a^^\u)) + a^^\u), where the noise component a^p{u) is the solution to 



Eq. (2.39 ) with the mechanical signal source term A{uj, uj') set to zero, while the signal 
component a^^\u) is the part of the solution to Eq. (2.39) that depends linearly on 



A{u,uj'). Thus, from Eq. (2.37) we can express the 'out' probe mode operator as 
follows: 



i^/27^e" 



+ 



^'r'^aP{uJ) + af{u) 



(2.46) 
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where the first square bracketed term gives the signal contribution to the detector 
response and the second square bracketed term gives the noise contribution. 

As 'in' states, we consider the mechanical oscillator bath to be in a thermal state 
at temperature T and the pump line to be in a coherent state centered about the 
pump frequency c^pifTD] 



\{a{uj)}), = exp dua{u) [a^+M - aJ^H] ) |0)„ (2.47) 
where |0)p is the vacuum state and 



aM^-M/^^-^. (2.48) 



The coherent state coordinate a{(jj) is parametrized such that the expectation value 
of the right-propagating 'in' current J'°(x, t) with respect to this coherent state has 
amplitude Jq, where 



r'^(x,t) = / duj./^ [e''^(^/''*'-*)aJ^(w) - e-*^(^/^''-*)a;^+(w)] , (2.49) 



with Vp = 1/ LpCp the wave propagation velocity in the pump probe hne. 

With the pump probe line in a coherent state, we assume that for large drive 



currents Eq. (2.39) can be approximately solved using a semiclassical, 'mean field' ap- 
proximation, where the quantum fluctuation Sa'^\u) in a^p{u)) = {a'^\u))) +Sa'^\u)) 
is kept to first order only. However, the nonlinear Duffing and transmission line- 
mechanical oscillator interaction terms can give rise to a bistability in the transmission 
line oscillator dynamics and one must be careful when interpreting the results from 
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the mean field approximation when operating close to a bifurcation point; large fluc- 
tuations can occur in the oscillator amplitude as it jumps between the two metastable 
amplitudes, which are not accounted for in the mean field approximation. (See Refs. 
[71] and [72] for respective analyses of the classical and quantum oscillator fluctua- 
tion dynamics near a bifurcation point). This issue will be further discussed in the 
following sections. 

The solutions to the signal (w) and noise a^'' terms parallel closely our previous 
calculations, which omitted the Duffing nonlinearity: [27] the Duffing {D) term in 



Eq. (2.39) has a very similar form to the transmission line oscillator coupling {B) 
term, both involving a|,aj operator combinations. We therefore relegate the solution 
details to the appendix, presenting only the essential results in this section. 

The solution to {a^^\uj)) is sharply peaked about the pump frequency oOp for large 
Tm and so can be approximately expressed as a delta function: {a'^\u})) = x^i^—^p)- 



Substituting this expression into Eq. (A. 3), we obtain for the amplitude x- 



X = c+ [2B{up, 0) + D{LOp)] X \xf , (2.50) 



with 



e='^^"7. (2.51) 
where Au = Up — Ux is the detuning of the pump frequency Up from the transmission 



line resonance frequency oot- Using the expressions for B{uip, 0) and D{ujp), Eq. (2.50) 
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can be written as 



{UJT -UJ,- ^7pt)x + ^^Xlxr = e'^^^^2TiIlZ^^pTl{fu^p). (2.52) 



where the effective Duffing couphng is defined as 



1^ = K,- ^^^K^m- (2.53) 



6m 



Notice that the interaction between the transmission hne and mechanical oscillator 
induces an additional Duffing nonlinearity (the second term involving in /C) in 



the transmission line mode amplitude effective equations of motion (2.52). However, 
in contrast with the Kd nonlinearity, which can be tuned to have either sign, the 
former mechanically-induced nonlinearity is always negative and thus has a "spring- 
softening" effect on the transmission line mode. Interestingly, by choosing an appro- 
priate compensating "spring hardening" Kd > 0, the effective Duffing constant /C can 
in principle be completely suppressed so that the next non-vanishing higher order 
nonlinearity would govern the mode amplitude dynamics. 

Once we have the solution for {a'^\uj)), the solutions for the quantum signal 
a'^\uj) and quantum noise 6a'^\u) are obtained from Eqs. (A. 2) and (A. 4), respec- 
tively. These solutions can be expressed as follows: 

aip\uj) = ai{ui)A{ui, u — ujp) + a2{io)A{uj — 2Auj, oo — ujp) (2.54) 



and 



5aP{u) = (3i{u)6C{u) + f32{uj)6C+{2ujp - u), (2.55) 
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where the ai{(jj) and Piiu) functions are defined in Eqs. (A. 8), (A. 9), (A. 12), and 



(A.13). 



Substituting Eqs. (2.54) and (2.55) into the expression (2.46) for a°"*(co') and then 



evaluating the signal component of the detector response (2.38), we obtain 



signal 



IoKTmUJT\'^ IpT 



7pT 



^ -[2niu-Up) + l] + 



c c \uj — Up ~ Au + i'jpT 

2'Jbm 



-\ln{yip — w) + 1] 



7pT 



•m JuJs~Sw/2 



271 \up{u-Up + AujY + 7p2. 



ai{(jj) ^ a2{(jj) f Lu - ujp + Auj + i'^pT 



sin [{u - uj')Tm/2] 

C \ijJ — Up — Aoj + ijpT J {u — Uj')Tm/2 
^Fext(c^ - UJp)F*^^{u' - Up) 



u -Up + Aw)2 + 7^ 

^ Fext C^p - u)F^^^[Up - U 

[up-u- umr + iL 



(2.56) 



where n{u) = (^^''^^/^bT _ ^^-i jg ^j^g thermal average occupation number for bath 
mode u. In the limit of small drive current amplitude Jq — > 0, we have ai{u)/c — 1, 
a2{u)lc — )■ 0, and we see that the signal spectrum comprises two Lorentzian peaks 
centered at Up^Um- The Up + Um peak corresponds to phase preserving detection, in 
the sense that a°"* gives the amplified a™ signal, while the Up — u^ peak corresponds 
to phase conjugating detection, with amplifying the a^^^ signal. [73] Increasing the 
drive current amplitude causes the peaks to shift, and the peak widths relative to their 
height to change, signifying renormalization of the mechanical oscillator frequency and 
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damping rate. The noise component of the detector response is 



1 ^ 



noise Zp J^^_g^/2 277 



-hw- 



pT 



|/32Mr-Re[/3i(u;)] 



UJ - ujp + Au;)2 + 7^^ 

[oj — ujp + Acj) 



IpT 



+Z. 



2 2ti'' 



Im[/3i(a;)] 
(2.57) 



where the integral term involving the fii{uj) functions includes the back reaction noise 
on the mechanical oscillator and the term involving Zp describes the probe line zero- 
point fluctuations added at the output. 



In Sec. 2.5 we will numerically evaluate Eqs. (2.56) and (2.57) and in particular 



compare the detector noise with the minimum noise bound discussed by Caves [271 [73] 
that follows from the Heisenberg uncertainty principle for the detector: 



mm— noise 



^ 2 27r \ ^pT J I'pT + Aco' 



X 



X 



ujs-Suj/2 27r yup {u - Up + Auy + -f^^ ^ 

'2'Jbm '^Ibm 



ai[uj) ^ a2[uj) f Lu - ujp + Auj + i'jpT 



UJ — Up — Auj + i'jpx 

(2.58 



2.4 Bistability Conditions 



We have seen [Hamiltonian (2.26)] that the current-dependent SQUID effective in- 
ductance gives rise to a transmission line Duffing type nonlinearity with strength Kd- 
Furthermore, there is a nonlinear coupling with strength i^rm between the transmis- 
sion line and mechanical oscillator. These two nonlinearities correspond respectively 
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to the cubic terms proportional to Kd and K"^^ in the mean transmission hne coor- 



dinate amphtude x equation (2.50). For sufficiently large drive current amplitude /q 



and/or coupling strengths KTm-, Kd-, the cubic term x \x\ term in Eq. (2.52) becomes 
appreciable, resulting in three real solutions over a certain pump frequency range Wp. 
This parameter regime defines the bistable region of the detector phase space (the 
intermediate amplitude solution is unstable and cannot be realized in practice). In 
the following, we determine the conditions on the parameters for the bistable region 
employing the analysis of Ref . [58] . 

We first express the transmission line mode coordinate in terms of its phase and 
amplitude: 

,(0) 



{af{t)) = Me-^('"^*+^"), 



X 



27iMe 



-i<Pm 



(2.59) 



where the amplitude M is a positive real constant and recall x is defined through the 
relation {a^p{uj)) = x^{^ ~ ^p)- Equation (2.52) then becomes 



(2.60) 



where 



2huJr, 



(2.61) 



Multiplying both sides of Eq. (2.60) by their complex conjugates and substituting 



E = M^, we obtain the following third-order polynomial in E: 



E^ + 



-E' + 



pT 



E 



(2.62) 
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The bifurcation line in current drive and detuning parameter space that dehneates 



between the single solution and bistable solution regions occurs where the susceptibil- 
ity dE/dujp diverges. If we further impose the condition that the transition between 
the two regions is continuous, i.e. d'^Up/d'^E = 0, we obtain the bistability onset 



critical point. From Eq. (2.62), these two requirements can be written 



SJC^E^ + A{uT - ujp)lCE + {UT - i^pf + YpT = 0, 
QK?E + A{ut - ujp)K, = 0. 



(2.63) 
(2.64) 



Solving these equations simultaneously for E and Au = Up — ut yields the following 
bistability onset critical values: 



E, 



Aubi = VS'jpT 



/CI 



Substituting these critical values into Eq. (2.62) gives 



?)\/?)UJT\fC\ 



Finally, using Eq. (2.61) we obtain the driving current critical amplitude: 



2fiuin 



'i\/?)UjT\KL\Zp 



(2.65) 



(2.66) 



(2.67) 



Note, the requirement that we operate below the Josephson critical current, Jq < /c, 
gives a lower limit on the value of |/C| for which our system can approach the bistability 
onset. The boundary of the bistable region that is given by the diverging susceptibility 
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equation (2.63) can be expressed in units of the bistability onset critical current // 



and detuning value Auu using Eqs. (2.67) and (2.65) to obtain [71] 



1 + 3 ' 



Im 2 V Aw, 



Aw 



± 



V Aw 



3/2' 



1/2 



(2.68) 



where the ± roots give the upper and lower boundaries of the bistable region, re- 



spectively (see Fig. 2.2). As mentioned in the preceding section, care must be taken 



when applying our semiclassical, mean field approximations to the detector signal and 
noise response when approaching closely the bifurcation boundary lines. Fluctuation- 
induced jumps between the small and large amplitude solutions of the transmission 
line mode can occur that are not accounted for in the mean field approximation. Nev- 
ertheless, in the next two sections we shall in some instances evaluate the detector 
response close the boundaries of the bistability region. For example, we shall see that 
significant improvements in cooling can be achieved provided a way is found to keep 
the transmission line mode on the low amplitude solution branch when operating in 



the bistable region. 



2.5 Displacement Detection 

Assuming that ^ 7^5-, i.e., the unrenormalized mechanical oscillator ampli- 
tude damping rate is much smaller than the transmission line oscillator amplitude 
damping rate, then the detector spectral noise and response in the mechanical signal 
bandwidth is approximately white over a large range of drive current and detuning 

31 



Chapter 2: Displacement detection and cooling of a nanomechanical resonator using 
a nonlinear microwave cavity detector 




Figure 2.2: Bistable region (shaded) of the cavity-oscillator system for negative, spring 
softening Duffing nonlinearity. The drive current and detuning are expressed in units 
of the bistability onset critical values lu and lAwtil- The labelled straight line traces 
correspond to detection {d) and cooling (c) current drive-detuning parameter exam- 
ples considered in Sees |2.5| and 2^ The arrows give the direction in which the drive 
current is varied in order to enter the bistable region on the small amplitude branch. 
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parameter space. The mechanical signal and noise response spectra are therefore 



approximately Lorentzian and Eqs. (2.56) and (2.57) can be parametrized as 



= Ga 



signal 

2'-fi,m[2n{R^u^) + 1] 



,S^I2 27r {u-UJpT RujUJmy + {R^lbr, 



(2.69) 



and 



{[Sn\tUs = tUp±R^tu^,SLo)Y) 



" 2mR,,,uj„ 



27back[2n±,k + l] 



271 ((jJ — =F Ru^n 



{R^lbr, 



out / 



Up ± RujUm, Su)] ) 



added noise 



(2.70) 



where G± is the phase preserving (conjugating) gain (in W-m~^), n{RujUJm) is the 
mechanical oscillator's external bath occupation number at its renormalized frequency 
RujOOm, R'yibm IS the reuormalized (i.e., net) mechanical oscillator damping rate, and 
the detector back reaction noise on the oscillator is effectively that of a thermal bath 
with damping rate 7back = IbmiR-y ~ 1) and thermal average occupation number 
^back- Note, from here on we do not consider an external classical force driving 
the mechanical oscillator; the focus is on displacement detection rather than force 



detection. The added noise term in Eq. (2.70) comprises output noise that is not due 



to the action of the detector on the mechanical oscillator; the added noise is present 
even when there is no coupling to the mechanical oscillator, i.e., when KTm = 0. In the 
absence of the transmission line Duffing nonlinearity, the added noise simply consists 
of the probe line zero-point fluctuations hwsSu/{4:7iZp). However, with the Duffing 
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nonlinearity present, the added noise will be in excess of the probe line zero-point 



fluctuations. 

The convenient Lorentzian parametrization approximations of the mechanical sig- 



nal (2.69) and noise response spectra (2.70) that provide the above-described effective 



thermal description of the back reaction noise will break down as one approaches arbi- 
trarily closely the jump points at the ends of the small or large amplitude transmission 
line oscillator solution branches occuring at the boundaries of the bistable region in- 



dicated in Fig. 2.2 This is a consequence of the diverging damping (i.e., ring-down) 



time of transmission line mode.[58j Thus, when numerically solving (2.56) and (2.57) 



to extract the effective thermal properties of the detector back reaction, it is impor- 
tant to always check the accuracy of the Lorentzian spectrum approximation. 

For sufficiently large gain (i.e., large current drive amplitude), we can neglect the 
added noise contribution and we have for the noise-to-signal response ratio when the 
mechanical oscillator external bath is at absolute zero [i.e., n{R^ujm) = 0]: 



signal 



7back 



(2.71) 



On the other hand, in the large gain limit the Caves noise lower bound (2.58) gives 



a noise-to-signal ratio of one. For large gain, we typically have |2?T,^g^^j^ + 1| ^ 1 and 
thus to approach the Caves bound necessarily requires |7back| ^ 7fem-[Zl] 

As an example, we numerically solve for the signal and noise contributions of the 



detector response, Eqs. (2.56) and (2.57) respectively, as well as the Caves lower bound 
on the quantum noise (2.58). We consider Duffing nonlinearities Kd = —3.4 x 10~^ 
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and Kd = (i.e., no nonlinearity) . The integrated signal and noise bandwidth is taken 
to be 600 = 2Ry'jhm- The corresponding parameter values are: probe line impedance 
Zp = 50 Ohms, transmission line mode angular frequency u)t/(2tt) = 5 x 10^ s^^, 
transmission line mode quality factor Qt = ujtI{'2.1pt) = 300, mechanical frequency 
C(;m/(2vr) = 4 x 10^ s~^, mechanical quality factor = ^m/{'^lbm) = lO'^, oscillator 
mass m = 10~^^ kg, Josephson junction critical current Ic = 4.5 x 10~® A, junction 
capacitance Cj = 10~^^ F, external flux bias $ext = 0.442 $0) aiid external field in 
the vicinity of the mechanical resonator i?cxt = 0.05 T. These values give a zero- 
point uncertainty A^p = 1.45 x 10^^'^ m and transmission line-oscillator coupling 
Krr. = 1.1 X 10^5. 
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Figure 2.3: Detector noise versus signal response at Aw = for harmonic [Kd = 0) 
transmission line. Duffing [Kd < 0) transmission line (di), and Caves' bound (black- 
dashed). Noise for the nonlinear transmission line is also evaluated for blue detunings: 
Au = +0.2\AuJbi\ (^2), and +0.4|Ac(;bj| (d^). The labeled curves correspond to the 
traces in Fig. 2.2[ The dashed, colored lines give Caves' bound for the corresponding 
detuning values. 
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The advantage of using a spring softening nonlinearity, Kd < 0, is clearly evident 



in Fig. 2.3, where we plot the noise versus response signal under increasing current 
drive for a transmission line both with and without Duffing term driven on resonance, 
Au = Up — ut = 0. We also plot the response of the nonlinear transmission line 
for several positively detuned values, Ao; = ujp — ujt > 0. Termination of the curves 
indicates the signal value at which the damping renormalization = 0, beyond which 
the derived solutions become unphysical due to the net mechanical damping rate 
becoming negative and hence the motion unstable about the original fixed point. Note 
that the same criterion, namely i?^ > 0, is employed throughout the paper in order to 
ensure stability of the system. Again, the semiclassical, mean field approximation is 
expected to break down in the vicinity of termination points, where large fiuctuations 



in the mechanical oscillator amplitude occur. In Fig. 2.3, we see that with positive 
detuning, we can further approach the Caves bound. However, this is at the expense 
of reduced gain; depending on one's point of view, large renormalizations of the 
mechanical oscillator damping rate (and frequency) due to detector back action may 
or may not be allowed in detector displacement sensitivity figures of merit, affecting 
the maximum achievable gain as one approaches more closely the Caves bound. 



The trends displayed in Fig. 2.3 can be partly explained by invoking Fig. 2.4, which 
indicates qualitatively the force on the mechanical oscillator due to the microwave 
transmission line 'ponderomotive radiation pressure' force, both with vanishing and 
with nonzero Duffing nonlinearity and also for 'red' and 'blue' pump frequency de- 
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tunings. The work done on the mechanical oscillator by the radiation pressure force 
during one period of motion, due to the delayed transmission line resonator response, 
is given by the area enclosed within the hysteresis loop[75l [76] and can be related to 
the steady-state back action damping rate through 



7back 



Wl 
r' 



(2.72) 



where W is the work done on the mechanical oscillator, E is the average oscillator 

energy and r is the period of motion. When frequency pulling is taken into account. 
Linear Regime Near Bistable Point 
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(a) Mechanical Oscillator Displacement (y) (b) 

Figure 2.4: Cartoon indicating the 'radiation pressure' force exerted on the mechanical 
oscillator by the transmission line mode during one cycle of mechanical motion: (a) 
the harmonic transmission line mode approximation; (b) approaching the onset of 
bistability. The work done on the oscillator is proportional to the area swept out 
during each cycle, considerably exaggerated here for clarity. Positive mechanical 
damping (red detuning) results on the positive slope side of the curves. Negative 
mechanical damping (blue detuning) results on the negative slope side. A spring 
softening nonlinearity can result in improved cooling for red detuning and improved 
signal-to-noise amplification for blue detuning. 



the usual notions of red-detuned (Aw < 0) or blue-detuned {Au > 0) hold only in 
the weak drive limit. We will assume red (blue)-detuned to correspond to drive and 
detuning values Au where the net work done on the oscillator is negative (positive) 
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as seen in Fig. |2.4[ For a harmonic transmission line and for low drive powers, the 



frequency pulling effects can be ignored, since the effective Duffing coupling Eq. (2.53) 
is proportional to the square of the transmission line-mechanical oscillator coupling 
Ktiu) which contributes only weakly for the considered parameter values. Conversely, 
the Duffing term causes frequency pulling even at low input power and can signif- 



icantly alter the slope of the response curve. From Eq. (2.72), the decreased slope 



on the blue detuned side leads to a decrease in the damping rate magnitude which. 



through Eq. (2.71), leads as demonstrated above to a closer approach to the Caves' 
limit. As mentioned above, benefits in lower noise-to-signal resulting from further 
detuning deep into the blue region are offset by diminished achievable signal gain 
levels. 



Tuning the sign of the Duffing coupling /C (2.53) to be positive, so that we have 
a hardening spring, results in an increased back action damping rate for blue detun- 
ing, and hence a corresponding decrease in signal to noise relative to the harmonic 
transmission line resonator detector case. 



2.6 Cooling 



Referring to the parametrizations (2.69) and (2.70) of the signal and noise com- 



ponents of the detector response, we define the mechanical oscillator's net occupation 
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number through the following equation: 

7„et(2n^et + 1) = 7bm [2n{R^uJm) + 1] + 7back(2n^^^k + 1)> (2-73) 

where the net damping rate is 7net = Ibm + 7back = R-^lbm- The oscillator's net 
occupation number is then 

2nnet + l = R-^ [2n{R^uj^) + 1] + (1 - R-^){2n^^^ + 1). (2.74) 

In order to cool a mechanical oscillator to its ground state using detector back action, 
we therefore require a large detector back action damping rate, equivalently large 
damping rate renormalization i?^ ^ 1, together with a small detector back action 
effective occupation number n^^^^ <^ 1. 



Referring to Fig. 2.4, operating closer to the bistability increases the negative 



work done per cycle on the oscillator by the cavity and hence increases the back 



action damping rate for given current drive. In Fig. 2.5, we plot the mechanical 



oscillator damping rate renormalization factor R^, using the same parameter values 



as in Sec. 2.5 (e.g.. Duffing coupling = —3.4 x 10~^), but with a larger yet 
still feasible mechanical quality factor Qm = 10^ (which we shall adopt throughout 
this section). We clearly see the enhanced damping as one approaches the onset of 



bistability given by hi (2.67) and AtUbi < (2.65). 



For the example parameter choices of Sec. 2.5, we have oJm/lpT ~ 0.5 and thus 
we are operating in the so-called bad cavity limit, where cooling close to the ground 
state (i.e., n^et ^ 1) is not possible. [271 IZZl [ZH] While it is not difficult to achieve 
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Detuning (lAcoJ) 



Figure 2.5: Mechanical oscillator damping renormalization factor for detunings 
both above and below the bistable detuning Auju- The amplification region corre- 
sponds to negative back action damping, i.e., < 1. 



the good cavity limit Um > IpT simply by realizing sufficiently large quality fac- 
tor superconducting microwave resonators, together with high frequency mechanical 



resonators, [30] it is nevertheless worthwhile to address how nonlinearities can im- 
prove on the cooling limits in the bad-cavity case. With the fundamental motivation 
to demonstrate macroscopic quantum behavior, the anticipated trend is to work with 
increasingly massive and hence lower frequency oscillators, making it progressively 
more difficult to achieve the good cavity limit. 



In Fig. 2.6, we plot the dependence of detector's noise effective back action oc- 



cupation number nback on microwave drive current amplitude at the detuning bias 
Au = - 



UJ. 



m + 7pT' where |Aa;| < |Aa;{,j|. This is the optimum detuning in the har- 



monic, transmission line oscillator approximation, i.e., when nonlinear effects are 
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ignored. The noise effective occupation number is indicated for both a nonzero 
{K(i = —3.4 X 10^^) as well as zero {Kii = 0) Duffing nonlinearity transmission 
line. We also show for comparison the effective back action occupation number when 
the frequency pulling effects of both the ponderomotive coupling and Duffing 
coupling are neglected. The latter case is obtained by dropping the nonlinear 



microwave mode amplitude term in the mean field equation (2.50). The sharp rise in 
occupation number and associated sharp drop in damping renormalization at larger 
current drives is a consequence of crossing over into the amplification region due to 
negative frequency pulling of the cavity response relative to the fixed detuning. The 
decrease in occupation number as / ^ is accompanied by weak back action damp- 
ing, which prevents cooling the mechanical oscillator to such occupation numbers. 
Note that at smaller current drives the damping renormalization in the presence of a 
Duffing nonlinearity peaks above the corresponding damping renormalization without 
the Duffing nonlinearity. This damping enhancement can be qualitatively explained 



with the aid of Fig. 2.4 In the presence of the nonlinearity then, improved cooling 
can be achieved for smaller current drives. 

According to the above discussion, any improvements in mechanical oscillator 
cooling are due solely to enhancements in the detector's back action damping rate 



for given drive; as can be seen from Fig. 2.6, the absolute minimum attainable de 



tector effective occupation number is the same both in the presence and absence of 
the transmission line resonator Duffing nonlinearity. While the effects of enhanced 
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Figure 2.6: (a) Detector noise effective back action occupation number versus current 



drive when red-detuned at Aw = — y + 7pr, |Au;| < |Ac<;f,j|, with a Duffing non- 

hnearity (sohd hne), without a Duffing nonhnearity (dashed hne), and without both 
Duffing and ponderomotive nonhnearities (dotted hne). (b) OsciUator couphng renor- 
mahzation factor for the corresponding back-action occupation number curves. 



These plots are obtained for the straight line trace labeled Ci in Fig. 2.2 



back action damping may be beneficial in situations where one is facing constraints 
on the maximum achievable drive power. [50] it would nevertheless be more signifi- 
cant if reductions in detector effective occupation number could similarly be achieved 
through nonlinear effects. To see how this might be possible, we consider detunings 
corresponding to the pump frequency being to the left and away from the cavity res- 
onance, i.e., |A(jj| > |Aw;,i|, A(jj < 0. For such detunings, the mechanical oscillator 
'sees' a transmission line resonator effective quality factor that is determined by the 



steeper slope on the left side of the response curve (see Fig. 2.4). As we drive the 



transmission line resonator towards the lower bistable boundary (see Fig. 2.2), the 
slope of the response curve increases sharply and mimics a resonator with larger qual- 
ity factor, effectively getting closer to the good cavity limit and hence resulting in a 
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lower detector occupation number. [271 EH EH] Continuing to drive the transmission 
line resonator into the bistable region, and assuming that the resonator can be main- 
tained on the low amplitude, red-detuned solution branch, [79] the detector effective 
occupation number further decreases while the back action damping rate on the me- 



chanical resonator increases (as explained by Fig. 2.4). Eventually, the transmission 



line resonator becomes unstable at the upper bistable boundary indicated in Fig. 2.2 



and the oscillator jumps to the larger amplitude, blue-detuned solution (see Fig. 2.7). 




-3-2-10 1 2 3 

Detuning (|Aaii,j|) 



Figure 2.7: Transmission line resonator response curve for Qt = 300 restricted to the 
small amplitude solution branch. The example drive currents from bottom-to-top: 
I /hi = 0.8, 0.95, 1.15, and 1.3. The jump between small (red-detuned) and large 
(blue-detuned) amplitude solutions is indicated by the dotted lines. 



In Fig. 2.8, we plot the dependence of the detector effective occupation number 



'T-back on current drive for an example detuning value of Aw = — 2^^;^ -|- = 
l.'i^uju- Driving the nonlinear transmission line resonator towards the upper bound- 



ary of the bistable region (see Fig. 2.4) produces a sharp decrease in detector oc- 
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Figure 2.8: (a) Detector noise effective occupation number versus current drive when 
red detuned at Acu = 1.3 Aw^j, corresponding to straight hne trace c'^ in Fig. 2.2 The 
Duffing nonhnear transmission hne resonator occupation number (sohd hne) rapidly 
decreases as the resonator is driven towards the upper bistable boundary, assuming 
the resonator can be maintained on the small amplitude metastable stable solution 
branch. In contrast, a harmonic transmission line resonator (dashed line) or a cavity 
with neither Duffing nor ponderomotive nonlinearities in its mean field microwave 
mode equations (dotted line) shows no such decrease in the occupation number, (b) 
Mechanical oscillator damping renormalization factor for the same fixed detuning and 
drive current range. The dashed vertical lines indicate the boundaries of the bistable 
region for the given transmission line resonator parameters. 
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cupation number, and an occupation number value of (2nback + 1) ~ 0.55 can be 
obtained, well below that achievable when ignoring frequency pulling effects. The 
harmonic cavity shows no such decrease in occupation number, indicating the quali- 
tatively different quantum dynamical dependencies on and and the necessity 
of the former. We can quantify the effect of frequency pulling by comparing with 
a harmonic transmission line resonator with a quality factor value chosen so as 
to give the same detector effective occupation number. For the occupation number 
value (2nback + 1) ~ 0.55, we have Qt? ~ 600, corresponding to Wm/Tpf = 0.95, 
and therefore the mechanical oscillator behaves as if it is coupled to a cavity with 
double the quality factor. This translates into lower net mechanical temperatures 



as shown in Fig. 2.9, where we give the net oscillator occupation number rinet (2.74) 
for various external bath temperatures. The combination of nonlinearly-enhanced 
coupling R^jbm and enhanced transmission line effective quality factor can be seen to 
significantly affect cooling of the mechanical motion, even for relatively large external 
temperatures. 



In the numerical solutions to Eqs. (2.56) and (2.57), the Lorentzian parametriza- 



tions (2.69) and (2.70) were found to give good approximations even when the upper 
bistable boundary is approached quite closely. This is a consequence of the wide sepa- 
ration in the relaxation rates that determine the line widths of the harmonic transmis- 
sion line resonator and unrenormalized mechanical oscillator modes, i.e., 7?,^ ^ 7pT- 
The upper bistable boundary has to be approached pretty closely in order for the 
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Figure 2.9: (a) Net mechanical occupation number at Aw = — yw^ + 7pT for a 

harmonic transmission line resonator. External bath temperatures: T = 1 (solid line), 
10 (dashed line), 50 (dotted line) and 100 (dot-dashed line) mK. (b) Dependence 
of the net mechanical oscillator occupation number on current drive for a Duffing 
transmission line with detuning Atj = 1.3 Acjf,j. The bistable region boundaries are 
indicated by the dashed vertical lines. 



nonlinear transmission line resonator ring-down time to exceed the renormalized me- 
chanical oscillator damping time, resulting in the breakdown of the effective thermal 



description of the detector back reaction. The deviation of Eqs. (2.56) and (2.57) 



from the assumed Lorentzian response near the upper bistable boundary gives a gen- 
eral criteria for measuring how close this limit can be approached. In all of the plots 
shown in this section, the Lorentzian approximation is a good one over the resolvable 
scale of the plots. The actual minimum temperature that can be achieved depends on 
the upper drive threshold where the Lorentzian approximation breaks down, as well 
as on the ability to keep the transmission line resonator on the small amplitude solu- 
tion branch; the latter condition becomes progressively more difficult to satisfy as the 
upper boundary is approached, owing to the increasing probability of noise-induced 
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jumps to the large amplitude branch. 

A Duffing transmission line resonator nonlinearity can also produce cooling gains 



in the good cavity limit. In Fig. 2.10 we consider a transmission line resonator 
with Qt = 1000, giving Um/'^pT = 1-6, and compare the nonlinear transmission 
line resonator with the harmonic resonator approximation at optimal harmonic de- 
tuning. Again, by detuning to twice the optimal harmonic resonator value, Au = 



— 2yc(;^ + 7py ^ 2.2Aubi, we see that the effective back action occupation num- 
ber decreases, while the back action damping increases as the system is driven to- 
wards the upper boundary of the bistable region. Driving a Duffing transmission line 
resonator at twice the optimal harmonic detuning can yield a detector occupation 
number (2nback + 1) ~ 0.06 just below the upper boundary of the bistable region, 
which is equivalent to an effective harmonic resonator quality factor of ~ 1400 
or Um/lpT — 2-2. In comparison, the minimum effective detector occupation number 



ignoring nonlinear effects is 2nback + 1 = 0.13. In Fig. 2.11 we plot the net mechan- 
ical occupation number for the good cavity transmission line resonator both in the 
presence and absence of the Duffing nonlinearity. Again, we see the strong cooling 
effects provided by frequency pulling of the cavity response. As discussed above, the 
minimum achievable net occupation number will depend on the threshold drive for 
which the Lorentzian approximation breaks down, as well as on the ability to lock the 
transmission line resonator onto the small amplitude solution branch in the bistable 
region. 
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Figure 2.10: (a) Detector noise effective occupation number versus current drive when 



red-detuned at Aw = — yi^m + TpT' l^'^l I^i^mI, for a Duffing nonlinear (solid 
line), harmonic (dashed line) transmission line, and with the effects of frequency 
pulling due to both ponderomotive and Duffing nonlinearities neglected (dotted line). 
The vertical dashed lines give the bistable region boundaries for the Duffing and 
harmonic transmission line resonators. This plot is obtained for the straight line trace 
labeled C2 in Fig. 2.2 (b) Oscillator coupling renormalization factor for optimal 
harmonic detuning, (c) Detector occupation number detuned at twice the harmonic 
optimum, Ao; = 2.2Au;f,j and locked to the lower stable amplitude solution. This 
plot is obtained for the straight line trace labeled in Fig. 2.2 (d) Corresponding 
back-action damping rate when driven to the upper bistable boundary. 
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Figure 2.11: (a) Net mechanical occupation number at Aw = - 

harmonic transmission hne resonator. External bath temperatures: T = 1 (solid 
line), 10 (dashed line), 50 (dotted hne) and 100 (dot-dashed line) mK. (b) Dependence 
of the net mechanical oscillator occupation number on current drive for a Duffing 
transmission line with detuning Au = 2.2Aubi. 



2.7 Conclusions 



We have provided a quantum analysis of a nonlinear microwave amplifier for dis- 
placement detection and cooling of a mechanical oscillator. The amplifier comprises 
a microwave stripline resonator with embedded dc SQUID. The SQUID gives rise to 
an effective, Duffing-type nonlinearity in the fundamental microwave mode equations, 
as well as a ponderomotive-type coupling between the microwave and fundamental 
mechanical modes. It was found that a spring-softening Duffing nonlinearity enables 
a closer approach to the standard quantum limit for position detection as expressed 
by the Caves bound, as well as cooling closer to the mechanical mode ground state. 
These findings can be qualitatively explained by considering the effects of frequency 
pulling in the response curve of the transmission line resonator 'ponderomotive force' 
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acting on the mechanical oscillator (see Fig. 2.4). With blue detuning, the decrease 
in damping allows for a closer approach to the quantum limit with large amplifier 
gain. Conversely, red detuning towards the bistable point of the force response curve 
increases the back action damping, improving the thermal contact to the detector 
'cold load'. Furthermore, effectively increasing the cavity quality factor due to the 
nonlinearity mimics the so-called good cavity limit in the harmonic case, allowing 
cooling closer to the ground state. 

The present investigation has by no means exhaustively searched the large param- 
eter space of the transmission line resonator-embedded SQUID-mechanical resonator 
system for establishing the optimal displacement detection sensitivity and cooling 
parameters. Rather, our intention has been to point out general trends, using specific 
parameter values as illustrative examples. It may be that other choices of parameters 
(e.g., using a mechanical resonator with a smaller quality factor) lead to a closer 
approach to the standard quantum limit, or to cooling closer to the ground state. 

The semiclassical, mean field methods employed in the present work do not take 
into account classical or quantum noise-induced jumps between the small and large 
amplitude metastable solutions that become more likely as the bistability region 
boundaries are approached. Unless ways can be found to keep the transmission line 
resonator locked onto the smaller amplitude solution branch, the predicted effects of 
nonlinearity-induced cooling will be less substantial, as it will be necessary to operate 
deeper in the bistability region to avoid jumps. The driven microwave mode ampli- 
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tude dynamics in the vicinity of the bistable region boundaries is still a relatively 
unexplored area that requires more sophisticated theoretical techniques in order to 
elucidate the fluctuations between the small and large amplitude metastable solution 
branches. [HI EQl EH [721 [821 [831 [H] This will be the subject of a future investigation. 
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Chapter 3 



Analogue Hawking radiation in a 
dc-SQUID array transmission line 



3.1 Motivation 

In the previous chapter we saw how the embedding of a dc-SQUID into a mi- 
crowave cavity results in a nonlinear Duffing oscillator with the strength of the non- 



linearity governed by the amplitude of an externally applied magnetic flux, Eq. (2.28). 
One can therefore modulate the dynamical properties of the microwave resonator by 
applying a time-dependent flux through the SQUID loop. In this chapter we wish 
to extend the SQUID-embedded cavity model to allow for control of the dynamics 
in both time and space. To this end, we now consider an effective 1-1-1 dimensional 
system involving a superconducting coplanar waveguide with centerline formed from 
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an array of dc-SQUID^ In this configuration a dc-SQUID approximated as a lumped 
inductor forms an LC oscillator together with the geometric capacitance of the trans- 
mission line ground planes. Therefore this setup is essentially an array of coupled 
oscillators each with a nonlinear fiux-dependent frequency. Quite remarkably, the 
propagation of free-field^ in curved space-times is very similar to a set of coupled 
oscillators with time- dependent frequencies [85| 186] suggesting we should be able to 
create a 1+1 dimensional analogue of photons propagating in an effective curved ge- 
ometry in this device. In this chapter we will show that this is indeed the case by 
considering perhaps the most important example of quantum field dynamics in curved 
spacetime: the Hawking effect [57]. 

3.2 Introduction 

The possibility of observing Hawking radiation [87j in a condensed matter system 

was first suggested by Unruh who uncovered the analogy between sound waves in a 

fluid and a scalar field in curved space-time [88] . In particular, the fluid equations of 

motion can formally be expressed in terms of an effective metric matching that of a 

gravitating spherical, non-rotating massive body in Painleve-Gullstrand coordinates 

1+1 dimensional system is one that can be described using a field theory comprising a single 
time and space dimension. 

free field interacts only with the classical gravitational background. 
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[89] 

ds"^ = -[cl- v{rf] dt^ + 2v{r)drdt + dr^ + r^dQ^, (3.1) 

where Cs is the speed of sound and v{r) is the spatially varying velocity of the fluid. 
For a sound wave excitation in the fluid, with velocity Cg, the horizon occurs where 
f^(r) = and the excitation is incapable of surmounting the fluid flow. Since Un- 
ruh's original proposal, Hawking radiation analogues have been proposed using Bose- 
Einstein condensates , liquid Helium , electromagnetic transmission lines [U^ , 
and fiber-optic setups [93]. Estimated Hawking temperatures in these systems vary 
from a few nano-Kelvin to lO^K respectively, far above temperatures predicted for as- 
tronomical black holes and thus usher in the possibility of experimental observation. 
Additionally, the understanding of the physics associated with laboratory system 
analogues may provide clues as to resolving unanswered questions associated with 
Hawking's original calculation such as the trans-Planckian problem [M] . 

In this chapter, we propose using a metamaterial formed from an array of direct- 
current superconducting quantum interference devices (dc-SQUID's). Modulation of 
the propagation velocity, necessary for the formation of an horizon, is accomplished 
through application of an external flux bias through the SQUID loops as indicated 



in Fig. 3.1 1. Under appropriate conditions, this configuration provides the supercon- 
ducting realization of Ref. [92], with the benefit of available fabrication methods. 
Indeed, arrays of SQUID's with parameters near those required to observe the Hawk- 
ing effect have already been constructed [95l|60]. Furthermore, as a quantum device, 
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Figure 3.1: a) Layout of the dc-SQUID transmission line. We assume each SQUID 
element is formed from identical tunnel junctions with critical current Jc and capac- 
itance Cj. b) Effective lumped circuit model valid for frequencies below the plasma 
frequency and negligible SQUID self-inductance. 



the SQUID array goes beyond the capabilities of previously proposed systems, allow- 
ing the possibility to probe the effect on Hawking radiation of quantum fluctuations 
in the space-time metric. Thus, in principle, this setup enables the exploration of 
analogue quantum gravitational effects. 



3.3 Model 

We consider a coplanar transmission line composed of a centerline conductor 



formed by a long, ^ 1, series array of dc-SQUID's indicated in Fig. |3.1^ . For 
simplicity, we assume that all Josephson junctions (J J) have identical critical cur- 
rent Ic and capacitance Cj values. For an individual dc-SQUID, with 0i and 02 
representing the gauge invariant phases across the JJ's, the equations of motion for 
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7± = =t </'2) /2 take the form 

' + -% + cos(7.) sin(7_) + ^ = (3.2) 



with plasma frequency cjp = (27rJc/Cj$o)^^^; characteristic frequency = '2t^IcRn/^o, 
and normahzed self-inductance Pl = 27rL/c/$o- The parallel, normal current resis- 
tance of the junction is denoted Rn, while $o = h/2e is the flux quantum and $ext is 
the external flux through the SQUID loop. If Pl then the SQUID dynamics can 
be approximated by a J J with a fiux-tunable critical current, = 2/c cos(7r$cxt/'!&o); 
the dynamics of which can be written 

^^ + sm(7+) = -, (3.3) 

where we have dropped the damping term, assuming the temperature is well below 
the superconducting critical temperature, and where the effective plasma frequency is 
given by Lo'p = ^2-kPJ{2Cj^q). We will assume the validity of this approximation and 
consider a flux-tunable array of Josephson Junctions (JJA). If we additionally restrict 
ourselves to frequencies well below the plasma frequency and currents below the 
critical current, then a JJ behaves as a passive, flux and current dependent inductance 
given by 

L„(J„, <|.^,J = ^ arcsin {IJI^) (3.4) 



for the nth JJ in the array. The equivalent circuit is given in Fig. 3.1| 3 where we have 



labeled the length and capacitance to ground of each J J by a and Co, respectively. 
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Using Kirchoff 's laws, we can write the discrete equations of motion as 



Vn+l ~ Vn — — ; In+1 ~ In — ~C*0 — (^-S) 



From (3.4), we see that by controlhng the external flux bias, or by creating a varying 
current in the transmission line, we are able to modify the inductance and thus 
propagation velocity inside the transmission line. Here, we focus on using the flux 
degree of freedom as our tunable parameter. Creating a space-time varying current 
pulse, as in Ref. ^3|, can also be accomplished in our device. However, our simplifled 
model does not admit the correct dispersion relation to support the required stable 
nonlinear solitonic localized pulses in the parameter region of interest. Charge solitons 
can however be produced in the high impedance regime of our device [96] . 

3.4 Effective Geometry and Hawking Temperature 

By deflning potentials An such that J„ = —CodAn/ dt and Vn = A^ — An-\ [92j, 



the equations of motion (3.5) can be combined to yield the discretized wave equation. 



^LnCo^An = An+i - 2An + A^-l. (3.6) 
dt dt 

For wavelengths much longer than the dimensions of a single SQUID the dispersion 
relation becomes to lowest order in k: 

u;\k) = -^sm'('^] ^c'e, (3.7) 
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where we have defined the velocity of propagation as c = a/ y/LCo, which in practice 
is well below the vacuum speed of light Cq. In this limit, the wave equation approaches 



the continuum 



did 92 
dt (? dt dx^ 



A = 0. 



(3. 



By ignoring higher-order terms in Eq. (3.7), we effectively remove the discreteness 



of the array which, along with dispersion from JJ inertia terms, can play the role 
of Planck scale physics in our system [071 ESI MM ■ For parameter values considered 
below, the relevant short distance scale is c/up (> a). Requiring the propagation 
speed to vary in both space and time. 



— )■ c^{x — ut), 



(3.9) 



with fixed velocity u set by an external fiux bias pulse. Fig. 3.2, the wave equation in 



flux pulse 



SQUID array 




bias line 



Figure 3.2: Layout of the dc-SQUID array transmission line used in generating ana- 
logue Hawking radiation. The fiux pulse creates a space-time varying magnetic field 
that modulates the speed of light in the SQUID array leading to an effective horizon 
in the comoving frame where u"^ = c^(x). 
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the comoving frame becomes 



d 



d 



dt ^ dx I \ dt dx 



d 



d 



a. n 



dx'^ 



A = 0, 



(3.10) 



where x and t now label the comoving coordinates. This wave equation can be re- 



expressed in terms of an effective space-time metric, 



eflf 



(? — 



-U 



-u 



(3.11) 



7 



Comparing this metric with Eq. (3.1), we see that our system contains a horizon 



located wherever = c^{x). In Fig. 3.3 we plot the effect of a step-like hyperbolic 



tangent flux bias pulse, similar to that shown in Fig. 3.2, with amplitude $ext = 0.2$o 



on a JJA with inductances given by Eq. (3.4), where we have kept only the lowest term 



in the Ic/Ic expansion. Additionally, since $cxt can only increase the inductance, the 




axis of propagation (x) 



Figure 3.3: Effect of a step-like flux pulse on the propagation velocity of a JJA as 
seen in the comoving frame. The pulse velocity was chosen to be m = 0.95c ($ext = 0). 
The black hole horizon occurs where c{x) = u. Arrow indicates the only permissible 
direction of travel across the horizon. 
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flux-bias pulse velocity u must be below the unbiased transmission line propagation 
velocity c in order to establish a horizon. We do not consider Gaussian or similar 
pulse shapes as they generate both black hole and white hole horizons [98] which 
complicates interpretation of the emission process. 

So far, we have focused on demonstrating a classical effective background geom- 
etry with an event horizon. The next step is to quantize small perturbations in the 
potential field A about this background. The correct commutation relations between 
quantum field operators are required for conversion of vacuum fluctuations into pho- 
tons [99]. These relations have been verified in the systems to which ours is analogous 
[52] . The resulting Hawking temperature is determined by the gradient of the JJA 
velocity at the horizon 

' . (3.12) 



Th 



2tx kh 



dx 



The radiated power in the comoving frame coincides with the optimal rate for single- 
channel bosonic heat flow in one-dimension [92| I100[ 1101] 



Eq. (3.13) is universal for bosons since the channel-dependent group velocity and 
density of states cancel each other in one dimension [100] . For a detector at the 
end of the transmission line, the radiation emitted by an incoming bias pulse will be 



doppler shifted yielding higher power compared to Eq. (3.13). However, the rate of 
emitted photons remains approximately unchanged. 
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For a single effective JJ, the magnitude of quantum fluctuations in the phase vari- 
able 7+ depends on both the ratio of Josephson energy, Ej = $o-^c/27i", to charging 
energy, Ec = e^/ACj, as well as on the impedance of the junction's electromagnetic 
environment. These energy scales give a representation of the phase-charge uncer- 
tainty relation A'yAQ > e, and relate the amplitude of quantum fluctuations between 
these variables [9]. When Ej/Ec > 1 and the impedance seen by the junction is less 
than the resistance quantum, Rq = h/Ae^ ~ 6.45kn, the phase operator behaves as a 
semiclassical quantity, i.e. the 7+ quantum fluctuations are small with respect to its 
average, and the JJ is in the superconducting state, allowing for a lumped inductor 
approximation. In the majority of experimental configurations, a single JJ is con- 
nected to probe leads with impedance ~ 50fi and as such is in the low-impedance 
regime Z/Rq <^ 1. In contrast, a JJA has an environment that comprises not only 
the input and output ports, but also all the other JJ's in the array. In this case, we 
can define an effective impedance as seen by a single junction to he Zj = Ze + Za 
where Ze is the environmental impedance of the leads and Za is the array impedance 
that, for frequencies below the plasma frequency, can be written as |102j 



where the last equality explicitly shows the dependence on the external flux bias and 
single junction parameters. Thus, even for a small energy ratio Ec/ Ej, the lumped 



3.5 Model Validity 




(3.14) 
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inductor model applies only when Za/Rq ^ 1- In Fig. 3.4 we show the dependence 
of array impedance on the external bias for fixed critical current Jc = 2 fiA 
and a range of experimentally valid capacitances to ground. As $ext — > '^'o/2, high 




Figure 3.4: Ratio of array impedance Za to the resistance quantum Rq as a function 
of the external flux bias for a critical current Ic = "2 /zA and example ground capacitor 
values: Co= 10-i6(solid), 5 x lO-^^(dashed), lO-i^(dash-dot), and 5 x lO-iSF(dotted). 
High and low impedance regions are defined above and below Za/Rq = 1 respectively. 



impedance causes large phase fluctuations, indicating a breakdown of our semiclassical 
description; the array undergoes a quantum phase transition from superconducting to 
insulating Coulomb blockade behavior [103J. Note, the small JJ parameter variability 



in actual arrays [951 iGOl ri03j will prevent the divergence in Fig. 3.4, as well as cause 
some transmission line scattering in the low impedance superconducting state. 

The dependence of the Josephson energy Ej on external flux $ext as described 
above allows for the systematic introduction of quantum fluctuations in our model. 
With the phase variable governing the circuit inductance, these fluctuations manifest 
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themselves in the effective metric (3.11) through the propagation velocity c. As the 
amplitude of fluctuations increases, the metric becomes a quantum dynamical variable 
which must be included in the description of the Hawking process. Thus, consequences 
of back-reaction from the Hawking process as well as quantum dynamical space-time 
can be probed by this configuration. Both processes, not included in the original 
Hawking derivation, represent analogue quantum gravitational effects present in our 
system |1U4] . 

3.6 Experimental Realization 



A possible realization of the JJA is shown in Fig. |3.5[ which consists of the JJA 
transmission line as well as an additional conducting line producing the space-time 
varying external flux bias $cxt- To provide a space-time changing velocity, the JJA is 
modulated by generating current pulses in the bias-line, the propagation velocity of 
which are assumed to be slightly below that of the unbiased JJA. The required bias 
pulse velocities u can be achieved by similarly employing individual JJ's in series as 
the bias line. Additionally, a dc-external flux can be used to fine-tune the transmission 
line velocity closer to that of the bias-line, eliminating the need for large amplitude 
current bias pulses. Unavoidable current pulse dispersion in the bias-line, resulting 
in a decrease in Hawking temperature, can be minimized with appropriate choice of 
pulse shape and transmission line length. 

Unambiguous verification of the Hawking process will require frequency-tunable, 
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tunable phase-qubit 
detector 




Figure 3.5: Possible transmission line and detector realization. Current pulses in 
the bias-line provide external flux necessary to modify the SQUID array propagation 
velocity. A phase-qubit at the end of the JJA functions as the photon detector. 

single-shot photon detection at the end of the JJA opposite to that of the bias pulse 
origin. Although not presently available, microwave single-photon detectors based on 
superconducting qubits are under active investigation [IS], I1U5] . We will assume a 
phase-qubit as our model detector [T8|. By repeatedly injecting current pulses down 
the bias-line, the predicted blackbody spectrum associated with the Hawking process 
can be probed by tuning the qubit resonant frequency. Correlations across the horizon 
between the emitted photon pairs can be established through coincidence detection. 
We emphasize the essential need for correlation information in order to establish that 
a photon is produced by the Hawking effect rather than some other ambient emission 
process, or spuriously generated via capacitive coupling to the bias-line. Unwanted 
directional coupling can be minimized with proper engineering of the transmission 
line. 

To estimate the Hawking temperature we will assume parameters similar to those 
of Ref. [95], with SQUID's composed of tunnel junctions with /c = 2 /iA and an 
upper bound achievable plasma frequency cjp = 27r x 10^^ Hz. The capacitance 



64 



Chapter 3: Analogue Hawking radiation in a dc-SQUID array transmission line 



to ground is assumed to be Co = 5 x 10 F (dashed line in Fig. 3.4). Using a 



SQUID length a = 0.25 fim gives an unbiased transmission line velocity c ~ Cq/IOO. 



Equation (3.12) gives the temperature as determined by the rate at which the JJA 
transmission line velocity varies that, in our case, is limited by the plasma frequency 
ojp. Assuming the maximum rate is an order of magnitude below ($cxt = 0) /2tt, 
then the Hawking temperature is ~ 120 mK; identical to a black hole with a mass of 
10^^ kg, or equivalently a Schwarzschild radius of 1.5 mm. This temperature can be 
a factor of ten larger than the ambient temperature set by a dilution refrigerator and 



therefore should be visible above the background thermal spectrum. Using Eq. (3.13) 



and the sample pulse in Fig. (3.3 ) gives an initial Hawking temperature 120 mK, which 



decreases ~ 10% every 1000 JJA elements due to bias-line dispersion. Applying the 



power expression (3.13) yields an average emission rate of one photon per pulse for 
~ 4800 squid's. Of course, the transmission line can be made considerably shorter 
at the expense of an increase in number of pulse repetitions in order to accumulate 
sufficient photon counts to verify the Hawking radiation. The parameters and pulse 
shapes chosen here illustrate feasibility of this setup, but do not represent the only 
available configuration. These values can likely be improved upon and optimized in 
terms of both performance and fabrication of this proposal. 
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3.7 Conclusion 

We have demonstrated that an array of dc-SQUID's in a coplanar transmission 
hne, when biased by a space-time dependent flux, creates an effective space-time 
metric with a horizon. As a quantum device, the superconducting transmission hne 
aUows for the possibihty of observing not only the Hawking effect, but also the effects 
of quantum fluctuations in an analogue gravitational system. 
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Chapter 4 

The Trilinear Hamiltonian: 
Modehng Hawking radiation from 
a quantized source 

4.1 Motivation 

We have shown that a magnetic field-pulsed microwave transmission line compris- 
ing an array of superconducting quantum interference devices, or SQUID's, not only 
reproduces physics analogous to that of a radiating black hole, but does so in a sys- 
tem where the high-energy and quantum mechanical properties are well understood 
and can be directly manipulated in the laboratory. This may help in addressing the 
validity of several assumptions made by Hawking as to the ultra-high energy physics 
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where quantum gravity effects are important. Missing from the analysis presented in 
Chapter [3] is the energy loss due to the emission of analogue Hawking radiation, and 
therefore this setup does not address the process of black hole evaporation. Indeed, 
the number of SQUIDs needed to produce a single photon per flux pulse is already 
pushing the current limits on experimental feasibility. However, hidden in conven- 
tional derivations of black hole evaporation is a breakdown of the unitary evolution 
in quantum mechanics known as the Information Loss Paradox[106j. Such a severe 
violation of the laws of quantum mechanics has stimulated a vast amount of literature 
devoted to solving this dilemma. Given that the evaporation process requires under- 
standing the laws of physics at the Planck scale, any solution must wait for a complete 
formulation of quantum gravity and is therefore likely out of reach for the foreseeable 
future. Building on our previous work of condensed matter analogues for black hole 
dynamics, in this chapter we will examine a simple nonlinear quantum optics Hamil- 
tonian as a toy model for the emission of Hawking radiation by an evaporating black 
hole. Although the system presented here is highly simplified, we will show that it 
is capable of reproducing many of the dynamical features thought to be present in 
the original problem. In addition, this work suggests that Hawking radiation may be 
entangled with gravitational degrees of freedom leading to nonclassical states of the 
black hole. 
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4.2 Introduction 

In the 35 years since Hawking's seminal paper on the quantum emission of radi- 
ation from a black hole |87j. a large body of work has been devoted to solving the 
so called information loss problem. For a black hole of fixed mass M, this emis- 
sion process yields a black body spectrum with characteristic temperature Th = 
hc^ /Snk^GM , irrespective of the initial state of the matter from which the black hole 
is formed. The inability to reconstruct the initial, possibly pure state of the black 
hole from the total emitted radiation signals the apparent breakdown of unitary evo- 
lution and the S-matrix description of the Hawking process. Although Hawking's 
calculations have since been verified in a number of ways |107] 11081 11091 1110[ 1111] , this 
breakdown at the foundation of quantum mechanics suggests that our understanding 
of black hole dynamics is not yet complete. 

The information loss problem rests on two key assumptions made in the standard 
derivation of Hawking radiation: (1) the perfectly thermal (i.e. mixed) character of 
the outgoing radiation; (2) the validity of this emission process over the lifetime of the 
black hole. The notion of thermal spectrum considered here is not that of a blackbody 
frequency spectrum, but rather the quantum thermal probability distribution defined 
by the temperature Th of a single mode (single frequency) of Hawking radiation. 
The traditional picture of the Hawking process leaves no room for deviations from 
this thermal distribution and thus breaks the requirement of pure-state— J-pure-state 
evolution enforced by unitarity. By itself, this process need not lead to information 
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loss since the information content of the black hole may be stored in entanglement 
between particle pairs created on opposite sides of the horizon [112^ I113j . Results 
from many-body theory suggesting that entanglement across a boundary scales with 
the area of the boundary lends credence to this view |1 14] . However, with the second 
assumption, the black hole causes a loss of entanglement and thus information. 

With the expectation that information must be conserved, many suggestions 
for resolving the information loss problem have been put forward. Current pro- 
posed solutions include long-lived and stable Planck-scale remnants |1 151 fTT6] . baby 
universes 1 1 1 7[ [TT8] . and the possibility of information escaping as non-thermal Hawk- 
ing radiation |119[ ll2Uj . In all of these proposals, corrections to the Hawking pro- 
cess manifest themselves when the black hole has evaporated to a size near that 
of the Planck length, Ip = ^JhGJ(? ~ 10^^^ m, at which quantum gravitational 
effects, neglected in Hawking's original analysis, are expected to play a role. In 
considering quantum states of the gravitational field, there is the possibility of back- 
reaction and entanglement of the radiating matter degrees of freedom with those 
of gravity p. 2 11 I122j . Although it is natural to consider a quantized gravitational 
field for the Hawking process, the current lack of a full quantum mechanical de- 
scription of gravity severely limits progress in directly addressing this scenario. In 
fact, exactly which degrees of freedom, if any, should be quantized is still subject to 
debate P^H| I124j . In this paper we investigate a simple, zero- dimensional quantum 
optics model of Hawking radiation that mimics some of the essential physics present 
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in the original information loss problem. 

As a zero-dimensional model we consider the following trilinear Hamiltonian: 

H = hwaa^a + hoJbb^b + hUcC^c + ihx (^ab^c'^ — a^bcj , (4.1) 

consisting of three harmonic oscillator modes with the frequency relation, Ua = Uf,+u}c- 
We will designate the modes as pump (a), signal (6), and idler (c) respectively. This 
Hamiltonian describes several quantum optics processes including frequency conver- 
sion, Raman and Brillouin scattering, the interaction of two-level atoms with a single 
mode resonant EM field, and is the full quantum generalization of the paramet- 
ric amplifier [nSl [I2ZI [Ml [HHl [1301 [Ml [132] • The connection between black 
hole radiance and parametric amplification was appreciated shortly after Hawking's 
discovery |1 33 j . Both processes amplify vacuum fluctuations resulting in the produc- 
tion of correlated photon pairs. Tracing over one of the two subsystems (i.e., signal 
and idler) yields statistics that are identical to that of a thermal distribution |134l [135] . 
The energy source in the parametric amplifier is assumed to be a classical pump 
such as a laser or microwave generator with fixed amplitude driving a system with 
X^^^ nonlinearity, the first nonlinear susceptibility in a medium without inversion 
symmetry |136j . Viewed as a black hole model, the pump plays the role of black 
hole mass M while the signal and idler modes of the parametric amplifier represent 
the escaping and trapped Hawking photons, respectively, as depicted schematically 



in Fig. 4.1 The trilinear Hamiltonian (4.1) generalizes the parametric amplifier by 



quantizing the pump mode and allowing for energy loss to the signal and idler modes. 
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Figure 4.1: a.) Diagram depicting the dynamics of the trihnear Hamiltonian. Initial 
vacuum modes are not shown for clarity. The nonlinear interaction is generated 
by a system possessing a second-order susceptibility x*-^"*- b.)Equivalent dynamical 
elements involved in the Hawking process. 



The expectation value of the pump mode energy is analogous to the mass M of a 
quantum mechanical black hole. We will explore this model by establishing the con- 
ditions under which the signal mode spectrum of the trilinear Hamiltonian deviates 
from the predicted thermal spectrum of the conventional parametric amplifier. We 
shall see that the quantization of the pump mode degree of freedom results over time 
in entanglement with the signal and idler modes and dynamics that become markedly 
different from the parametric approximation. In particular, the signal mode devel- 
ops a strongly non-thermal spectrum that is dependent on the initial pump mode 
state. The corresponding entropy is reduced relative to that of a thermal (maximally 
mixed) state indicating the presence of information. These model system results lend 
support to the view that late-time Hawking radiation contains an increasing amount 
of information about the initial quantum state of the black hole and is composed of 
particles entangled with quantized gravitational states. 



The outline of this chapter is as follows: In Sec. 4.3 we review the derivation of the 
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amplification of vacuum fluctuations under the parametric approximation. Sec. |4.4 
considers the semi-classical approximation whereby the back-reaction from the quan- 
tized radiation onto the classical pump is accounted for and derives the self-consistent 



equations of motion for the signal and idler modes. In Sec. 4.5 we consider the full 



quantum dynamics Eq. (4.1 ) under the short time analytical approximation and com- 



pare this result along with the full quantum numerical solution to the parametric and 



semi-classical approximations from the previous sections. Sec. 4.6 investigates the 



role of entropy and entanglement in the production of non-thermal states. In Sec. 4.7 



we give a possible microwave cavity realization of the trilinear Hamiltonian. Finally, 



Sec. 4.8 concludes with a brief discussion of the results and consequences for black hole 



evaporation. Details of the numerical calculations are presented in Appendix B.2 



4.3 Parametric Amplifier and Hawking Emission 

In the following, we derive the well-known thermal spectrum of the signal mode 
under the parametric assumption of a fixed amplitude pump mode. Replacing the 



pump mode in Eq. (4.1) with a fixed amplitude drive A results in the interaction 
frame Hamiltonian 

Hi = ihxA (6+C+ - be) . (4.2) 
The Heisenberg equations of motion for the signal and idler mode operators are, 

^ = A^c(t)-: '-f = A^m*. (4.3) 
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These can be readily solved yielding the Bogoliubov transformations, 

6(r) = 6(0) cosh (v4r) + c(0)+ sinh (Ar) (4.4) 
c(r) = c(0) cosh (At) + 6(0)+ sinh (At) , 

where we have expressed the dynamics in dimensionless time r = xt- If the system 
starts with both signal and idler modes in the groimd state, |\1'(0)) = |0,0)^^, then 



Eq. (4.4) gives for the number operators A^^ and N^. 

iVb(r) = N^{t) = sinh^ (At) . (4.5) 

Additionally, we are interested in the probability distribution of the individual signal 
and idler subsystems in the number state basis. With the system initially in the 



ground state, the unitary evolution corresponding to Eq. (4.2) can be expressed as, 

|^(r)) = exp [At (6+c+ - be)] |0, 0)^^ . (4.6) 
Making use of the disentangling theorem |137j . 



exp 



[At {b+c-^ - be)] = ^rb+c+^-g{bn+c+c+i)^-rbc^ ^^ ^^ 



with r = tanh (y4r) and g = In cosh (Ar), where the last term in Eq. (4.7) vanishes 
for the ground state and the middle term reduces to e~^, we have 

e'""'" |0,0),, = I0'0)^c = Y.'^nh-{AT) |n,n),,, (4.8) 

n=0 ■ n=0 

with the evolving state vector |\l'(r)) given by 

|^(r)) =sech(Ar)Etanh"(Ar)|n,n)j,^. (4.9) 

n 
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Let us now focus on operators acting on a subsystem spanned by the states of either 
the signal or idler mode individually. Assuming we are interested in the signal mode 



only, tracing over the idler subsystem in Eq. (4.9) gives the operator expectation value 



(0(r))j, = sech^ (At) ^ tanh^" (At) {n\0\n)^ (4.10) 



for an arbitrary signal mode operator O. Comparing Eq. (4.10) with the spectrum of 
a thermal state defined by temperature T, 

n n 

indicates that the signal mode is in a thermal state provided we define the temperature 



as 



^^^^"2A;Bln[coth(Ar)]' ^^'^^^ 
where the time dependence is a consequence of the rapid increase in the occupation 



number of the modes, Eq.(4.5). As in the Hawking process [138j . the particle pairs 



generated by the parametric amplifier form an entangled two-mode squeezed state 



given by Eq. f4.9) |139j . The bipartite structure of this system allows calculating the 
entanglement between particle pairs via the Von-Neumann entropy 5*^, also referred 
to as entanglement entropy. For mode i = b,c, using the reduced density matrix pi 
we have 

5, = -Tr(panpi), (4.13) 
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where we drop the usual fee factor. For both the signal and idler modes, this entropy 
is given by 

which is the thermal entropy of a quantum harmonic oscillator with temperature 



defined by Eq. (4.12). Thus, as for the Hawking thermal radiation, we see that the 
temperature of the parametric oscillator is determined by the entropy generated from 
tracing over one of the two modes in a particle pair squeezed state. 



4.4 Semi- Classical Analysis: Backreaction on a Clas- 
sical Pump Mode 

We now go beyond the just-considered parametric amplifier model and incorporate 
backreaction effects due to the emission process using a semi-classical approximation 
where the pump mode is treated as a classical variable, while the signal and idler 
modes are quantized. In order to self consistently solve for the system evolution, we 
first work out the dynamics of the pump mode assuming it behaves as a classical 
variable affected by the expectation values of the signal and idler modes, after which 
we substitute the resulting c-number expressions for the pump operators into the 
Hamiltonian and calculate the signal and idler evolution. A similar procedure arises 
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for the semi-classical Einstein equations: 



G,. = ^(^\T,^\^), (4.15) 



c4 



where the mass M black hole spacetime geometry is modeled classically through 
the usual Einstein tensor G^^, interacting with the quantum matter/radiation fields 
via the expectation value of the stress-energy-momentum tensor operator T^^ with 
respect to a suitable incoming quantum field state |\&). Solving these semiclassical 
equations using an analogous procedure to that outlined above results in a nonlinear 
Schrodinger equation jl4U] and also a possible loss of spacetime stability |141j . In 
addition, issues such as non-causal dynamics |142] and possible inconsistencies with 
wavefunction collapse |143j occur. In contrast, the semi-classical model presented here 
can be readily solved by exploiting the symmetries present in the Hamiltonian. 



To begin, we consider Eq. (4.1) in the interaction frame 



Hi = ihx{ab+c^ -a^bc) (4.16) 



and obtain the following mode equations, 

da , db , dc , , 

dt = rft = ' dt= ^^-'^^ 

that lead to the evolution of the pump mode number operator, 

^ = -^ (a6+c+ + a+bc) . (4.18) 

It is easy to show the corresponding number operators for signal and idler modes 
are given by dN^/dt = dNc/dt = —dNa/dt. To proceed, we will use the following 
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Manley-Rowe constants of motion [144j . 



Mab = Na + Nh; Mac = Na + N,; Mbc = Nh-N, 



(4.19) 



expressing the underlying SU(2) and SU(1,1) symmetries in our modelfHHl 1146] . 



Differentiating Eqs. (4.17) again results in decoupled equations of motion containing 



only commuting operators. 



dt^ 
d^Nh 

df^ 



2x^ [3iV2 - Na (2M,b + 2Mac + 1) + M,,M,J 

-2x' [3iV2 - N, (4M,, - 2Af„, - 1) + M,, {M,, - M,, - 1)] (4.20) 

-2x' [3iV2 - {AMac - 2Mab - 1) + Mac (Mac ' M^b " 1) 



We now take the expectation value of the pump number operator equation (4.20) 



for Na and make the approximation (A^"^) ~ (Na) (Na), the validity of which can be 



checked using full quantum numerical simulations of Eq. (4.17) [see Sec. 4.5 . These 
approximations, along with the condition that both signal and idler mode start in 



the vacuum state, lead to the semiclassical evolution for the pump mode: 



iV„(r) = /3+ + [iV,(0) - dn 



iVa(O) - /3- 



(4.21) 



where dn{u, k) is the Jacobi elliptic function and 



/3: 



1 + 



2iV,(0) ± v/l + 12iV,(0)+4iV,(0)2 



(4.22) 



It is important to note that both Eqs. (4.21) and (4.22) are expressed in terms of 



the initial conditions of the pump mode only, a consequence of the relations (4.19) 
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and the fact that the signal and idler are initially in the vacuum (ground) state. 
Equations of motion for both signal and idler can then be obtained by substitution 



of the c-number expressions for the pump mode with amplitude given by Eq. (4.21) 



(4.23) 



where 0(t) is a slowly varying function of time, into Eq. (4.16): 



(4.24) 



The time evolution of this now bilinear Hamiltonian can be straightforwardly calcu- 
lated if we assume |130 



dt 



const, 



(4.25) 



which is the most general condition under which the Hamiltonian at different times 



commutes, 
relation 



Hi{t),Hi{f) 



0. Combining Eqs. (4.17) and (4.23) results in the phase 



,da da~^ d(h(t) 
a+— - a— = -2iNa- ^ ^ 



dt dt 



dt 



(4.26) 



which, expressing the left-hand side as Eq. (4.16) through the Heisenberg equations 



for the operators Eq. (4.17), yields a second condition for the pump phase: 



2iV,(t) 



d(f){t) 
dt 



HA. 



(4.27) 



As long as the initial state contains at least one mode in the vacuum state, we have 



the constant of motion {Hi) = 0. Thus, Eqs. (4.25) and (4.27) are only consistant if 
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(l){t) = 0. This leads to a bilinear Hamiltonian of the form, 

Hj = thxNa{ty/^ {b+c+ -be) . (4.28) 



Using a similar analysis to that employed in Sec. 4.3 gives 

db{T) 



dr 
dcjr) 
dr 



6(0) cosh {9{t)) + c+(0) sinh {9{r)) (4.29) 
c(0) cosh (0(r)) + 6+(0) sinh (^(r)) , 



where. 



Oir) = / ViVa (r')cir'. (4.30) 



Therefore, we see that each mode behaves similarly as in the parametric approxima- 
tion but with the number of particles dependent on the pump dynamics: 

Nbir) = N,{t) = sinh^ (^(r)) . (4.31) 

However, most importantly the modes are again in a thermal state with distribution 

(0(r)) = sech^ (^^/nJ^)t^ ^tanh^" (^^/lU^)r^ {n\0\n) . (4.32) 

n 

Thus, the spectrum of the signal mode (Hawking radiation) remains thermal when 
pump mode energy loss due to signal and idler particle creation is taken into account 
within the semi-classical approximation. 
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4.5 Full Quantum Description 

4.5.1 Short-time Approximation 

In this section we consider the full quantum dynamics of the trilinear Hamiltonian 



Eq. (4.1), (i.e. quantize the pump mode as well). Although the trilinear system does 
not allow for a general analytical solution, we can obtain approximate equations in 
the short time limit r = ^ 1- Using the condition Ub = oUc = uJa/2, appropriate 



for modeling particle generation by a black hole, we can rewrite Eq. (4.1) ignoring 
constant terms as 

H = Ho + Hi^t = {a+a + K,) + ihx {aK+ - a+K_) , (4.33) 



where 



K+ = b+c+, = be, K, = ^ {b+b + c+c + l) . (4.34) 



We now construct the Casimir operator. 



= Kl-]^ {K^K_ + K_K^) = i [Ml - l) , (4.35) 



where the second equality comes from definition (4.19). This operator obeys the 
eigenvalue equation 

K'^\^) = k{k-l)\^) , (4.36) 

where k = 1/2 {\Mbc\ + 1). We now look for simultaneous eigenstates of both 
and Kz'. \k,nc), where Uc is the number of idler particles. These states can also be 
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decomposed using the number state basis 

\k,nc) = \nc + 2/c - 1)^ 1^^)^ = \nb) In^) , (4.37) 

with Hfy representing the number of particles in the signal. Operating on this state 
with Kz gives 

Kz \k,nc) = {k + He) \k,nc) = | + + M |^^^ |^^^ ^ (4.38) 



valid only for initial conditions where both signal and idler modes are in pure states 
satisfying > ric- Including the pump mode state in the full state vector and 
expressing the idler mode population as a function of the pump amplitude 

\na)\k,Mac-na) , (4.39) 

one may switch to the interaction frame where the evolution depends only on Hi^^ 



(4.33) 



|r; k, Mac) = en-^+-+^-) |0; k, M„,) , (4.40) 
with the initial state with Mac = ^a(O) is given by 

\Q-k,Mac) = \Mac)\k,Q). (4.41) 



The short time approximation r = ^ 1 to Eq. (4.40) can be calculated using the 



Baker-Campbell-Hausdorff formula truncated to O (r^) 

^T(^aK+-a+K^) ^ ^TaK+ ^-Ta+ K- ^- ^[aK+,a+ K.] ^ (4.42) 
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Acting with the third exponential operator term on our initial state gives 



(4.43) 



indicating clearly that the time over which the approximation is valid t = r/x < 



^/ixV^^Mac) decreases with pump amplitude and coupling strength. Evaluating the 



BCH approximated expression (4.42) on the initial state (4.41), we obtain for the full 
evolution of the state: 



Ma 



|r; k, Mac) = , V /„ (fc, Mac) r" \Mac -n)\k, n) 

where NM^d'^) is the time-dependent normalization factor 



(4.44) 



NMaAr)=e-~'T'''-T{Mac + l,r-') 



(4.45) 



with r (a, b) the reduced gamma function and 



fn{k,Ma 



MaJT {2k + n) 
n\ {Mac-n)\T{2k) 



1/2 



(4.46) 



Our interest in the crossover from classical to quantum dynamics for the pump mode 



suggests that we use coherent states built from linear combinations of Eq. (4.44) for 
the initial pump state. To this end, we consider a general initial state 



1^(0)) = x; 



a.^M4=|.)|0)|0), 



s=0 



(4.47) 



with the pump mode in an as yet unspecified initial pure state with probability 

2 

Ps = \as\ of being in state s and both signal and idler modes in the vacuum state 
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{k = 1/2). Here we implicitly assume the probabilities are normalized and add to 



unity. The density matrix at some later time r resulting from (4.47) is given by 



,=0 .=0 .=0 ,=0 

which, after performing a partial trace, leads to the reduced density operators 

Pa (^) = E E E E ^sK. J^^ip^s, ^ I . - z) (r - J I , (4.49) 
where 5ij is the Kronecker delta, and 

Pnr) = j:±Ps§^^r^^\^^\ (4.50) 

s=0 i=0 '^^ ' 



for the pump and signal modes respectively. Fig. 4.2 compares the expectation values 



of the pump and signal modes for the parametric, semi-classical, and short-time quan- 



tum approximations, as well as the full numerical solution to Eqs. (4.17) where the 
pump mode is initially in a coherent state with amplitude (Ai"a(0)) = 9 corresponding 
to classical pump modes with amplitude A = 3, and the signal and idler modes are 
initially in their vacuum (ground) states. 

4.5.2 Mode Spectrum Dynamics Under the Short-time Ap- 
proximation 

The main limitation of the short-time approximation is the inability to account 
for backreaction effects resulting from the build up of quanta in the signal and idler 
modes. For particles produced in the Hawking process, however, the entangled pairs 

84 



Chapter 4-' The Trilinear Hamiltonian: Modeling Hawking radiation from a 
quantized source 




(a) 



0.4 0.6 O.f 
Time (t) 

(b) 



Figure 4.2: a.)Evolution of pump mode initially in coherent state (iVa(O)) = 9 for the 
parametric (red), semi-classical (blue), short-time quantum approximation (green), 
and full quantum numerical solution (black), as a function of dimensionless time 
r = xt- b.) Corresponding population of the signal mode. Note that the semi-classical 
analysis ceases to be valid when the pump mode becomes depleted. 



generated early in the evolution do not, to first-order, effect those created later from 
the black hole [106]. Furthermore, the emitted radiation does not build up in the vicin- 
ity of the black hole, but escapes to spatial infinity. Therefore, in this section we will 
suppose that the expressions derived above for the evolving pump and signal states. 



Eqs. (4.49) and (4.50) respectively, in fact provide a more relevant zero-dimensional 



model of a radiating black-hole when extrapolated beyond their original short-time 
domain of validity. The pump and signal modes contain in general a large number 
of coefficients for each number state and cannot be easily evaluated. However, in the 
long-time limit where the pump is nearly depleted, these equations can be consid- 



erably simplified by noting that a general element of Eq. (4.50) for fixed s is given 



85 



Chapter 4-' The Trilinear Hamiltonian: Modeling Hawking radiation from a 
quantized source 



by 



T{s + 1) 



u = s — I. 



(4.51) 



'Ns{t)' ''M!r2«r(s + l,r-2)' 
As such, only those elements where i = s remain nonzero, leading to the signal mode 

density matrix 

oo 

P* = E^^I^)(^I' (4-52) 

s=0 

which is in general a mixed state with number-state probabihty distribution iden- 
tical to that of the initial pump state. Therefore, by measuring the late-time signal or 
idler mode distribution we recover all but the phase information associated with the 



initial pure state of the pump mode. Additionally, Eq. (4.52) shows that the number 



of quanta in the signal mode is equal to the initial number of quanta in the pump 



mode, indicating that the pump (4.49) ends up in the vacuum (ground) state. 



The most important knowledge gained from Eq. (4.52) is that the signal mode 



spectrum will no longer be that of a thermal state, in contrast to the parametric 



and semi-classical approximations in Sec. 4.3 and Sec. 4.4, respectively. Focusing 



on coherent states, in Fig. 4.3 we give an example of the evolution of Eqs. (4.49) 



and (4.50) by plotting the probability distributions for both modes as a function 



of time r for a pump initially in a coherent state with amplitude (Ai"a(0)) = 9 and 
initial vacuum state for both signal and idler. Additionally, in Fig. |4.3) 3 we highlight 
what the thermal distribution would be at each time step by equating {Ni^ij)) with 
a Bose-Einstein distribution to extract an effective temperature: 



{Nk{T)) = [e^b/kMr) -I] \ 
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Figure 4.3: a) Spectrum of the pump mode at several time-steps when the pump is 
initially in a coherent state with (A^a(O)) = 9. b) Corresponding spectrum of the 
signal mode. Red lines indicate what the distribution would be at e ach ti me-step if 
in a thermal state corresponding to occupation number (A^f,(r)), Eq. (4.53). 



In Fig. 4.3 1 see we the evolution of the initial coherent state as it loses quanta to the 



signal and idler modes and evolves towards the ground state represented by a peak 
in the probability distribution at the n = number state value. The corresponding 



evolution of the signal mode in Fig. |4.3|3 starts with the vacuum state and progresses 



towards the state with probability distribution identical to that of the initial pump 



coherent state, Eq. (4.52). By comparison with the effective thermal state (4.53) 



we can see that the initial probability distribution for the signal mode is nearly that 
of a thermal state up until the pump mode has transferred nearly half of its initial 



energy corresponding to {Ni,{t)) = 4.5. As the evolution continues, Fig. 4.3 d shows 
the increasing deviation from the effective thermal description for the signal mode 



distribution as expected from Eq. (4.52) 



87 



Chapter 4-' The Trilinear Hamiltonian: Modeling Hawking radiation from a 
quantized source 



4.5.3 Non-thermal Spectra and Information 



We now quantify the deviations of the signal mode spectrum Eq. (4.50) from that 
of a thermal state using the fidelity [TTT] 



F,{t) = TT^Jp,{Ty/^a{T)p,{Ty/^, (4.54) 
where ^^(t) is the density matrix of the signal mode and cr(r) is a thermal density 



matrix with effective temperature determined by Eq. (4.53) using the occupation 
number {Nb (r)). The fidelity provides a measure of the distance between these two 
states in distribution space with range < F < 1, where unity indicates that the 



two density matrices are identical. In Fig. 4.4 1 we plot the fidelity of the signal mode 



assuming a pump mode initially in a coherent state starting with several occupation 



numbers ranging from one to nine. By comparison with Fig. 4.2, we can see that for 



the initial coherent state {Na{0)) = 9 the fidelity remains essentially unity, indicating 



that the signal mode density matrix Eq. (4.50) is nearly thermal, until the pump has 



transferred close to half of its initial energy into the signal and idler modes beyond 
which point there is a strong deviation from the thermal state as the signal mode 



asymptotically approaches the state given by Eq. (4.52). This is in agreement with 



the qualitative description presented in Fig. |4.3|3 and remains true for all the initial 



states considered in Fig. 4.4 1 



This deviation from a thermal distribution also indicates that the signal mode 
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0.0 0.2 0.4 0.6 0.8 
Time {r} 

(b) 



1.0 1.2 1.4 



Figure 4.4: a) Fidelities of the signal mode given a pump mode initially in a coherent 
state with occupation number (No) = 1 (blue), 3 (red), 6 (green), and 9 (black) 
under the short-time approximation, b) Information present in the signal spectrum 



calculated from Eq. (4.55). Equivalently colored vertical lines denote the time at 



which the effective subspace dimensions satisfy, d*^ = dl^. 



spectrum contains information defined as |119] 



h{T) = SI\t) - S,{t), 



(4.55) 



where (t) is the Von Neumann entropy of the signal mode in a thermal state with 



equal average occupation number, Eq. (4.14), and Sbir) is the actual entropy of the 



mode calculated using Eq. (4.13) and the reduced density matrix of the signal mode 



P{,(r). In Fig. 4b we plot the information contained in the signal mode for the initial 
pump coherent states considered in Fig. 4a. Given that the fidelity is nearly unity 
up until the pump mode transfers half of its original quanta to the signal and idler 
modes, it follows that the information content of the signal or idler mode is close to 
zero over this same time interval before increasing as the signal spectrum becomes 
identical to that of the initial coherent states. 
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In order to account for the dynamics of the signal mode information, we first 
consider a general bipartite pure state of a system with fixed total energy that is 
composed of subsystems A and B, each with finite Hilbert space dimensions dA and 
ds, respectively. It is known that subsystem B will be nearly thermal and thus con- 
tain approximately no information as long as dA ^ dB, with the information content 
of subsystem B becoming apparent only after dA^ ds |119t I148[ 1149] . Similar dynam- 



ics for the information content of the signal mode is shown in Fig. 4.4 However, the 



three harmonic oscillator modes from which a pure state of the trilinear Hamiltonian 



Eq. (4.1) is composed all contain an infinite set of states, preventing the direct appli- 
cation of dimensional arguments to our model. One can, however, define an effective 
subspace dimension for each mode p^9] 

^f(^) = Tr[aKr)] ' ' = ^^-^^^ 
determined by the purity of the effective thermal state cTj(r) with temperature given 



by Eq. (4.53). This definition is motived by the fact that the purity of a mixed state 
density matrix is proportional to the number of basis states over which the fractional 
population of the mixed state is nonzero. For a state with an energy (quanta) con- 
straint, using the thermal state cr(r) gives the minimal value for the purity, and as 



such Eq. (4.56) represents an effective maximum number of states constrained by the 
number of quanta in the mode at time r. Therefore, if we partition our initial pure 
state into bipartite subsystems da and composed of pump and combined signal 
and idler modes respectively, then the initial effective subspace dimensions satisfy 
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d'f S> dl^ = [df^Y , where the equahty is due to the symmetry between signal and 
idler modes. In Fig. 4.4 d we plot the times at which d'f^ = d'^ for each of the initial 



pump mode coherent states considered in Fig. 4.41. Just as for finite dimensional pure 
states, the information contained in the signal or idler subsystems remains nearly zero 
until after d'^^ir) ^ d^ij-), provided we define the effective subspace dimension using 



Eq. (4.56). Our results are in agreement with those of Page |119] where a similar 
argument was put forth for the evolution of information in the Hawking radiation 
from a finite dimensional black hole. 

4.5.4 Numerical Results for the Full Trilinear Evolution 

Unlike Hawking radiation, which is well-modeled using the quantum short-time 
approximation, backaction due to the build up of quanta in the signal and idler modes 



has a strong effect on the evolution of the trilinear Hamiltonian (4.16). In Fig. 4.2 we 
see that backreaction effects quickly lead to differing evolutions between the short- 
time and full quantum dynamics; backreaction prevents the full transfer of quanta 
from the pump mode [T50] . and results in the oscillation of quanta between pump and 
signal/idler modes. To account for these backaction terms we repeat the analysis in 



Sec. 4.5.3 by numerically calculating the full dynamics of Eq. (4.17). Fig. 4.5 shows 
the modifications to both the fidelity and information content of the signal mode 
when backaction is included in the dynamics, as well as the subspace dimension 



condition (r) = (i^^(r). As expected, the short-time dynamics in Fig. 4.5 are 
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Time ( 

(a) 



(b) 



Figure 4.5: a) Fidelities of the signal mode using numerical simulations of Eqs. (4.17) 
with a pump mode initially in a coherent state with occupation number (Na) = 1 
(blue), 3 (red), 6 (green), and 9 (black), b) Corresponding information present in the 
signal mode spectrum and times at which d'f^r) = dl^i^r) (dashed lines). 



nearly identical to those in Fig. 4.4 up until d'^{T) ^ dl^{T) which occurs later 
than in the short-time approximation as the backaction from particles in the signal 
mode begins to impede the transfer of quanta from the pump. The peaks in the 
information content of the signal mode correspond to the times at which backaction 
from the pump mode has completely stopped the transfer of energy between modes. 
With the majority of quanta now in the signal and idler modes, the flow of energy 
reverses directions as the signal and idler now drive the pump mode leading to the 



oscillations seen in Fig. 4^ As a model for black hole evaporation these oscillations 
represent the unphysical process of Hawking radiation flowing back into the black 
hole; the connection between the trilinear Hamiltonian and Hawking emission is valid 
only for the initial transfer of quanta from pump to signal/idler modes. 
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4.6 Tripartite Entanglement 

With the pump mode quantized, one can consider the entanglement between pump 



and signal/idler modes. In Sec. 4.3, it was shown that ignoring this entanglement 
(pump mode treated classically) and tracing over one mode of a two-mode squeezed 
state leads directly to the thermal characteristics of the remaining system. In the 
full quantum description the statistics of the signal mode is obtained by tracing over 
both idler and pump modes and as such the distribution of entanglement between 
modes of this now tripartite system is important for characterizing the spectrum of 
the signal mode. Additionally, entanglement with the pump mode may appreciably 
alter the state dynamics of the pump compared to the classical approximations in 



Sees. 4.3 and 4.4 Given that multipartite entanglement is not as well understood 
as for the bipartite case, we begin by again considering the system to be partitioned 
into two subsystems consisting of the pump mode and the combined signal-idler. 
This bipartite separation into pump and signal-idler subsystems allows us to use the 
mutual information |147j 

la-bc = Sa + She — Sabc (4-57) 

as a measure of total correlations, both classical and quantum, between subsystems |151j . 
Since our tripartite system state remains pure throughout its evolution, the total en- 
tropy Sahc = and the subsystem entropies satisfy Sa = Sbc- The mutual information 
is therefore twice the entropy of the pump mode subsystem, la-bc = 2S'a. Thus, the 
subsystem entropy of the pump is a direct measure of entanglement with the signal 
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and idler modes, which is not taken into account in either the parametric or semi- 



classical solutions. Of course the signal and idler subsystems are also entangled with 



each other as in Sees. 4^ and AA With both signal and pump modes in identical 
states, the signal-idler mutual information is given by 



'6-c 



Sb ~\~ Sc Sl)c — StSfo Sa 



(4.58) 



indicating the tradeoff between the entanglement of the signal/idler and that of the 



pump mode Eq. (4.57). In Fig. 4.61 we plot the mutual information Eqs. (4.57) and 



(4.58) for a pump initially in a coherent state of amplitude {Na{0)) = 9 for the para- 
metric, semi-classical and short-time approximations along with the full numerical 



solution obtained using Eqs. (4.17). We see the mutual information (pump entangle- 



ment) la-bc begins to become appreciable around the same time as the information 



content of the signal mode becomes apparent in Figs. K4b and 4.5 d indicating the 



increasing role of the pump mode in the dynamics. The numerical solution shows 
that this increase in pump entanglement reduces the signal-idler mutual information 
Ib-c with respect to the semi-classical and parametric approximations as expected 



from Eq. (4.58). In the full quantum dynamics the pump is never allowed to be fully 



depleted and is therefore always entangled with the signal and idler. However in the 
short-time approximation the late-time pump mode is depleted and approaches the 
ground state where Sa = 0. 

Finally, we show that the entanglement buildup with the pump mode given by 



' a~bc 



not only affects the signal and idler modes, but also results in amplitude de- 
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pendent squeezing of the pump. In Fig. |4.6{ 3 we plot the amphtude of the squeezing 
parameters [139] , 

b=4(A.Y|>-l. X, = "t^' ^;"'"> (4.59) 

for the pump mode, where (AX|.) is the variance of the quadratures. As can be 
seen, the entanglement with the signal and idler modes leaves the pump mode in a 
non-classical squeezed state over the time range of interest. 





0.2 0.4 0.6 0.8 1.0 
Time(r) 

(b) 



Figure 4.6: a) Signal-idler mutual information /f,_c for pump with {Na{0)) = 9 in 
the parametric (red), semi-classical (blue), quantum short-time (green), and full nu- 
merical solution (black). Dashed lines represent mutual information la-bc indicating 
entanglement with the pump mode, b) Quadrature squeezing parameters q± of the 
pump mode given the pump is initially in a coherent state with occupation number 
{Na) = 1 (blue), 3 (red), 6 (green), and 9 (black). Negative values indicating squeez- 
ing of the g_ quadrature. Dashed lines represent the corresponding g+ quadrature. 
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4.7 Circuit QED realization of the trilinear Hamil- 
tonian 



In attempting to fabricate a direct realization of the trilinear Hamiltonian Eq. (4.1 ) 



we run into difficulty when generating the specific nonlinear coupling between oscilla- 



tors. This coupling must not only obey the symmetries of Eqs. (4.19), but must also 
be tunable so that the signal and idler modes may be decoupled while the pump mode 
is driven to its initial energy state. Although in principle nothing prohibits us from 
creating such a coupling, at present we do not have a suitable model for this interac- 
tion. Instead we will make use of the connection between the trilinear Hamiltonian 
and the Dicke model [T25] in realizing our system as an array of superconducting 
qubits interacting with a single resonant mode of a microwave cavity 11531 1154] . 
In keeping with the theme of this thesis we will assume the use of fiux qubits formed 
by biasing a SQUID with an external magnetic fiux with amplitude near a half-fiux 



quantum, $ext = ^o/2. An illustration of this setup is given in Fig. 4.7 Of course, 
one may readily substitute one of the many other qubit varieties [12]. 

The Dicke Hamiltonian for a collection of qubits interacting with a single mode 
of a microwave cavity resonator is given as [125j 

H = hwcC^c + hoJcJz - ihx {c^ + c) (J+ - J-) , (4.60) 

where ojc is the frequency of the field mode, assumed to be resonant with the level 
splitting of the qubits, x is the coupling strength of the qubits to the cavity field, and 
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Figure 4.7: Cartoon of a typical coplanar transmission line resonator with flux qubits 
(SQUIDs) situated in between the centerline conductor and ground plane. In mod- 



eling the Dicke Hamiltonian, Eq. (4.60) the amplitude of the resonant cavity mode 



(red arrows) is assumed to be constant over the length of the qubit ensemble. 

{J±,Jz} are collective qubit operators satisfying the standard angular momentum 
commutation relations. 

[J+, J_] =2/iJ„ [J±,J,] = tJ± (4.61) 

It is implicit in this Hamiltonian that the only qubit interactions are with the cavity 
field and that the wavelength of this field is much larger than the size of an individual 
qubit; each atom experiences the same amplitude of the resonant cavity field. In 
the absence of any qubit-cavity interaction, the operators c(t), J-it) and c^(t), J+{t) 
rotate in the Heisenberg frame as exp {—iut) and exp {iut) respectively. Therefore, 
in making the rotating-wave-approximation (RWA) we can ignore terms containing 
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cJ- and c"*" J+ which gives the simphfied Hamiltonian 

H = hWcC^c + hoOcJz + ihx {cJ+ - c^J+) ■ (4.62) 

This Hamihonian can be rewritten using the Schwinger representation [155j where the 
angular momentum operators are replaced by products of bosonic operators 

J+ = a+6, J_ = a6+, Jz = ]^ (a+a - 6+6) , (4.63) 

with the boson operators a, a+ and 6, 6+ satisfying the usual boson commutation 
relations. These operators correspond to the annihilation and creation of excited and 



ground-level qubits respectively. The angular momentum operator J_|_, Eq. (4.63) 



now represents the switching of a qubit from the ground to excited state while the 



reverse process is governed by J_. Upon substitution of Eqs. (4.63) into Eq. (4.62) 
we obtain the Hamiltonian 

H = huatt'^a + hubb'^b + hUcC^c + ihx (a+fec + afe+c"*") (4.64) 

satisfying the frequency relation Ua = —ujb = which is equivalent to our trilinear 



Hamiltonian Eq. (4.1 ). The number operators A^i = a+a and N2 = b^b corresponding 
to the excited and ground state modes do not represent the actual mode population 
values but rather give the effective occupations obeying Ni + N2 = 2J < N, where J 
is the total angular momentum of the qubits. These effective values coincide with 
the actual level populations when all of the qubits are in either the ground or excited 
state, 2 J = N. 
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In constructing the trilinear Hamiltonian from a simplified Dicke model we have 
replaced the conceptually simple coupled cavity system with a much more abstract 
realization based on the effective oscillatory modes of the collective dynamics of a 
qubit ensemble. While at first this may appear disadvantageous, the connection 
to the Dicke model suggests that the dynamics of the information content in the 



signal and idler modes presented in Sees. 4.5.3 and 4.5.4 may be related to the well- 



known super-radiant phase transition of the Dicke Hamiltonian |156[ 1157^ 1158] . This 
connection will be addressed in a later work. 



4.8 Conclusion 

In this chapter we have investigated the effect of a dynamical quantized pump 
mode on the generation of quanta in a parametric amplifier. We have shown that 
a quantized pump mode leads naturally to a non-thermal spectrum for the signal 
and idler modes once the pump has released nearly half of its initial energy, such 
that the effective subspace dimensions of the pump and signal/idler mode systems 
approximately coincide. The departure from a thermal state indicates that the signal 
spectrum contains information that may be used to partially reconstruct the initial 
state of the pump mode. Once quantized, the pump mode becomes entangled with 
the signal and idler leading to non-classical squeezed states of the pump. As a simple, 
zero-dimensional model of the Hawking effect, the present findings lend support to 
the possibility for non-thermal emission from a quantum (as opposed to semiclassical) 

99 



Chapter 4-' The Trilinear Hamiltonian: Modeling Hawking radiation from a 
quantized source 

black hole; the emitted radiation is entangled with the quantized gravitational degrees 
of freedom and yields information concerning initial formation of the black hole. 
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Derivation of signal and noise 
expressions for SQUID detector 



In this appendix, we give the derivation of the signal a^r^\uj) and noise a')p\uj) 
terms of Eq. (2.46) in Chapter |2] Suppressing the signal dependent term A{u,u') in 



Eq. (2.39), we obtain the noise equation 



duj'B{uj,uj')a^\uj — Lo') 



-oo 
oo 



du" 



/oo roo 
du" / du' a^^^^ {u")a^\u')a''rp\u + u" — u') 
-oo J — oo 

+ C{u). 
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Keeping only terms to first order in A{uj,uj') or in Eq. (2.39), we obtain the signal 



equation 



X 



+ 



-oo 
oo 



du'a^^ {lo — uj')A{uj, u') + / duj' B{uj, c<;')a^'* [u — lo') 



duj" 



-oo 
oo 



du' B{u,u')aP {u — u') I du" a^''(a;")a^''^(t<;" — w') 



oo /"OO 

// 



+ D{uj) I duo" I duo' 

I (0)+/ ii\ (1)/ i\ (0)/ I // /\ 



(0)+/ //X (0)/ /X (1)/ , // /X 



(A.2) 



Decomposing aP{uj) = {aP{u)) + 6ap\uj) and expanding Eq. (A.l) to first order in 
the quantum noise fiuctuation 6aP{u), we obtain the following two equations: 



duj'B{uj,uj'){aP{uj -co')) 



'OO 

oo 



duo" 



/OO POO 
du" / du'{aP^iu")){aPiu')){aP{u + tu" - tu')) 
■oo J —oo 



(A.3) 



and 



du" 



du' B{u,u')5aP [u — u') 



+ 



■oo 
oo 



{aP{u")){ap-^{u" -u')) + {aP^{u")){aP{u" + u')) 



du'B{u,u'){aP{u -u')) / du" \6aP{u"){aP^{u" - u')) 
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+ 5aP{uj" + uj'){aP^{uj")) 

/OO P'QO |- 

duj" dtu' Sa^^^^{Lo"){a^^\tu')){a^^\LO + tu" -u')) 
OO J —OO 

+ 6C{lj). (A.4) 



Assuming {aip\u)) can be expressed approximately as a delta function, i.e., {a'i^'{u)) 



(0), 



xS{uj — Up), Eq. (A. 3) reduces to Eq. (2.50) for x- The semiclassical approximation to 
Eq. (A. 2) for a'>^\uj), with a^^(ci;) replaced by {a^p{uj)) = x^{^ ~ ^p)^ then becomes 



{1-2 Ixl' [B{uj, 0) + B{uj, uj-ujp) + D{uj)]} 

-X^ [2B{uj, uj -ujp)+ D{uj)] a^^^^{2ujp - u) = u -ujp). (A. 5) 

In order to invert and obtain a^j){uj), we require a second, linearly independent equa- 
tion that also depends on a^rp^(2ujp — oo). Such an equation is obtained by making 



the replacement to 2ujp — a; in Eq. (A. 5) and then taking the adjoint: 



{1 + 2 Ixl^ [B{uj - 2Aw, 0) + B{lo - 2Auj, to - cUp) + D{uj - 2Auj)]}a^^^^{2ujp - 

[2B{uj - 2Auj,uj- Up) + D{u - 2Au)] a^^\u) = -x*A{u - 2Au,u - Ujl)A.6) 



Now inverting, we obtain 



aip\u) = ai{u)A{u, u — Up) + a2{u)A{u — 2Au, u — Up), 



(A.7) 



where 



ai(cj) = V{u)-^ {1 + 2 Ixl^ [B{u - 2Au, 0) + B{u - 2Au, u - Up) + D{u - 2Au)] ] x, 

(A.8) 
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a2{uj) = -V{uj)-' [2B{u, u-ujj,) + D{uj)] \xf X 



(A.9) 



and the determinant is given by, 



V{uj) = {l-2\x\^[B{io,0) + B{io,io-ujp) + D{uj)]} 

X { 1 + 2 Ixl^ [B{uj - 2Auj, 0) + B{uj - 2 Aw, u - cUp) + D{uj - 2Auj)] } 

+ \xf[2B{uj,uj -Up) + D{u)][2B{u -2A(j , CO — (jjp) 

+ D(w-2Aw)]. (A.IO) 



A similar approach is used to obtain 5a^^\(jj) from Eq. (A.4), giving 



(A.ll) 



with 



= V{ujy^ {1 + 2 [B{u - 2Auj, 0) + B{u - 2Au,uj - cOp) + D{uj - 2Aw)] \xf} 

(A. 12) 

and 



/32(w) = V{u)-^ [2B{u, u-Up) + D{u)] x^- 



(A.13) 
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Included in this appendix are the source codes for the numerical simulations used 
in Chapters |2] and |4} All functions are coded for use in the MATLAB computing 
environment. In addition, the functions used in simulating the trilinear Hamiltonian, 
Sec. B.2 require the Quantum Optics Toolbox by Sze M. Tan |159] . 



B.l Signal, Noise, and Cooling Codes 

%Simulation codes for results obtained in 

%Phys. Rev. B 78, 104516 (2008) 

%Paul D. Nation, Dartmouth College, 2008 

%alpha 1 

function [out]=:alphal(w,IO , dw , Ktm , Kd , vars ) 
% alpha 1 in Eq.56 

out=det er (w, 10 , dw , Ktm, Kd, vars)." — l.*(l + 2.* (B(w— 2*dw, ,Ktm, vars ) + . . . 
B(w— 2*dw , w— (dw+vars ( 2 ) ) ,Ktni, vars ) + . . . 

D(w— 2*dw ,Kd, vars ) ) . * abs ( ChiNumeric (dw, 10 ,Ktni, Kd, vars)). "2).*... 

ChiNumeric (dw, 10 ,Ktm,Kd, vars ) ; 

%alpha 2 

function [out] = alpha2(w,I0 ,dw,Ktm,Kd, vars ) 
%alpha2 in Eq.56 

out=— deter (w, 10 , dw, Ktm, Kd, vars).' — l.*(2 . *B(w, w— (dw+vars ( 2 ) ) ,Ktm, vars ) + . . . 
D(w,Kd, vars )).*... 

abs ( ChiNumeric (dw, 10 ,Ktm, Kd, vars))"2* CliiNumeric (dw, 10 ,Ktm, Kd, vars ) ; 

%B 

function [out] ~ B(w, wprime , Ktm, vars ) 
%Eq.41 

out = ((vars (2) *Ktm/ sqrt(2*pi)) "2/2)./ (w— vars (2)+i*vars(4)).*... 

( 1 . / ( wprime— vars (3)+i*vars(5)) — l./( wprime+vars (3)+i*vars(5))); 

%hetal 

function [ out] = betal (w, 10 ,dw,Ktm,Kd, vars ) 
%betal in Eq.57 

out=det er (w, 10 , dw , Ktm, Kd, vars)." — l.*(l + 2.* (B(w— 2*dw, ,Ktm, vars ) + . . . 
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B(w— 2*dw , w— (dw+vars ( 2 ) ) ,Ktni, vars )+D(w— 2*dw , Kd, vars )).*.. . 

abs ( ChiNumeric (dw, 10 .Ktm.Kd, vars ) ) " 2 ) ; 

%heta2 

function [ out] = beta2 (w, 10 ,dw,Ktni,Kd, vars ) 
%beta2 in Eq.57 

out=deter (w, 10 ,dw,Ktni,Kd, vars).' — l.*(2. *B(w,w— (dw+vars ( 2 ) ) ,Ktm, vars ) + . . . 
D(w,Kd, vars ) ) . * ChiNumeric (dw, 10 ,Ktni,Kd, vars ) " 2 ; 

%G 

function [ out ] =C( 10 , dw, vars ) 
%Eq.42 

out=abs ( i * sqrt (2*pi)./(vars(4)— i . *dw ).*... 

sqrt (( 10 . " 2* vars (6)* vars (4) ) . / ( vars ( 1 ) . * (dw+vars ( 2) ) ) ) ) ; 

%ChiNumeric 

function [ out] = ChiNumeric (dw, 10 ,Ktm,Kd, vars ) 
%Upper or lower solution to Eq.50 in bistable region 
%determined by input vars—NLswitch 
if vars(9) ==0 

out=abs(C(IO ,dw, vars ) ) ; 

else 

p3 = l; 

p2 = -2.*dw./( vars (2)*K(Ktm,Kd, vars ) ) ; 

pl = (dw."2+vars (4) 2) ./( vars (2) " 2*K(Ktm,Kd, vars ) " 2) ; 

pO=-vars (4)/( vars (2) "2*K(Ktm,Kd, vars ) ' 2 ) . * . . . 

(10. "2* vars ( 6 ))/( vars ( 1 )* (dw+vars (2))); 
poly=ccll (1 , length (10 ) ) ; 
for j =l:length(IO) 

poly{l,j} = [p3,p2,pl,p0(l,j )]; 

end 

r=ccll (1 ,length(IO )); 
for j =l:length ( 10 ) 

r{l , j}=roots (poly{l , j }) ; 

end 

reals=zeros ( 1 , length ( 10 ) ) ; 
for j =l:length ( 10 ) 

reals (1 , j )=min( r { 1 , j } (imag( r { 1 , j }) = = 0)); 

end 

out=sqrt (2* pi ) * sqrt (reals ) ; 

end 

%D 

%Eq.43 

function [ out]=D(w,Kd, vars ) 

out = (vars (2 ) *Kd) / (2*pi).*l./ ( w— vars (2)+i*vars(4)); 

%deter 

function [out] = deter(w,IO ,dw,Ktni,Kd, vars ) 
%Determinant , Eq.AlO 

out = (l-2*abs ( ChiNumeric (dw, 10 ,Ktm,Kd, vars ) ) ' 2 . * (B(w, ,Ktm, vars ) + . . . 
B(w,w— (dw+vars (2) ) ,Ktm, vars)+D(w,Kd, vars ))).*... 

(l + 2*abs ( ChiNumeric (dw, 10 ,Ktm,Kd, vars ) ) ^ 2. * (B(w-2*dw,0 ,Ktm, vars ) + . . . 
B(w—2*dw,w— (dw+vars ( 2 ) ) ,Ktm, vars)+D(w— 2*dw , Kd, vars ))) + .. . 
abs ( ChiNumeric (dw, 10 ,Ktm,Kd, vars ) ) " 4 . * . . . 
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( 2 . *B(w, w— (dw+vars ( 2 ) ) ,Ktm, vars)+D(w, Kd, vars )).*.. . 

(2*B(w— 2*dw,w— (dw+vars (2) ) ,Ktm, vars)+D(w— 2*dw,Kd, vars ) ) ; 

%dwc 

function [ out]=dwc(Ktm,Kd, vars ) 
%Eq.65 

out=sqrt (3)* vars (4)*sign (K(Ktni, Kd, vars ) ) ; 

%G 

function [ out]=G(w, 10 ,dw,Ktm,Kd, vars ) 
%Gain function defined by Eq.56 

out=IO ."2./(2*pi).* (Ktm* vars (2)/ vars (4)) " 2* vars (4) "2/ ( vars (4) ~2+dw ' 2 ) . * . . . 
(w./ (dw+vars (2) ) . * vars (4) " 2 . / ( (w-(dw+vars (2) ) +dw) . " 2 + vars (4)"2)).*... 
abs ( alpha 1 (w, 10 , dw,Ktm,Kd, vars)./C(IO ,dw,vars) + ... 

alpha2 (w, ID ,dw,Ktm,Kd, vars ) . /C(IO ,dw, vars ) . * ( (w— (dw+vars (2))+dw+ . . . 

i*vars(4))./ (w— (dw+vars (2) ) — dw+i *vars (4)))). "2; 

%Ihi 

function [out] = Ibi (dw,Ktm,Kd, vars ) 
%Bistable current Eq.67 

out =2* vars (4)*sqrt((2*vars(l)* (dw+vars (2)))/(3*sqrt(3)*... 
vars (2)* abs (K(Ktni, Kd, vars))*vars(6))); 

%intminnoise 

function [ out] = intminnoise (w,dw, 10 ,Ktm,Kd, vars ) 
%integrand of min . noise (caves bound) in Eq.58 

out = 1/(2* pi ) . * (w/ (dw+vars (2)).*vars(4)"2./( (w- (dw+vars (2) ) +dw) . " 2 + . . . 

vars(4)"2)).* abs ( alpha 1 (w, 10 , dw , Ktm , Kd , vars )./C(IO ,dw, vars ) + . . . 

alpha2 (w, 10 ,dw,Ktm,Kd, vars ) . /C(IO ,dw, vars ) . * ( (w— (dw+vars (2))+dw+ . . . 

i*vars(4))./ (w— (dw+vars (2) ) — dw+i *vars(4))))."2.*... 

(2* vars (5) . / ( (w-(dw+vars (2)) - vars (3)) . "2 + vars (5) ' 2) - . . . 

2*vars(5)./((( dw+vars (2) ) — w— vars (3))."2+vars (5) ~2)); 

%intnois e 

function [out] = intnoisc (w, dw , 10 , Ktm, Kd, vars ) 
%integrand in Eq.57 

out=l/vars(6)*(vars(l).*w)./(2*pi)*2*vars(4)"2./( (w— (dw+vars (2) ) +dw) . " 2 + , . . 
vars (4) " 2) . * (abs (bet al (w, 10 ,dw,Ktm,Kd, vars ) ) . " 2 + . . . 
((w-(dw+vars (2) ) +dw) . " 2+ vars (4) " 2) . / ( (w-(dw+vars (2)) -dw) . " 2 + . . . 
vars (4) " 2) . * abs( beta2 (w, 10 ,dw,Ktm,Kd, vars ) ) . ~ 2 - . . . 

real(betal(w,IO, dw,Ktm,Kd, vars ) ) + (w— (dw+vars (2) ) +dw) ./ vars (4).*... 

imag (bctal (w, 10 , dw,Ktm, Kd, vars ) ) ) ; 

%K ■ 

function [out] = K(Ktm,Kd, vars ) 
%Effective nonlinearity , Eq.53 

out=Kd-(2*vars (2)* vars (3) *Ktm" 2)/ ( vars (3) "2+ vars (5) " 2) ; 

%lowbnd 

function [ out] = lowbnd (dw) 

%lower boundary of bistable region, Eq. 68 

out = (l/2).*dw. ~(3/2).*(l + 3.*dw."(-2)-(l-(l./dw)."(2))."(3/2))."(l/2); 

%minnoise 

function [ out] = minnoisG (ws , gamma, 10 ,dw,Ktni,Kd, vars ) 
%min . noise (caves bound), Eq.58 
int=@(x) intminnoise (x ,dw, 10 ,Ktni,Kd, vars ) ; 
out=abs (1/ vars (6)*(vars (1)* ws*2*ganima) / (4*pi) — ... 
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(IO*Ktm* vars (2)/ vars (4)) " 2* vars (4) " 2/ (dw~2+vars (4) " 2) * . . . 

quadgk ( int , ws— gamma, ws+gamma) ) ; 

%n 

function [ out]=n(w, vars ) 

%Bose— Einstein occupation number 

out = (exp ((vars (l).*w)./(vars(7)*vars(10)))— 1).' — 1; 

%n oise 

function [ out] = noise (ws , gamma, 10 , dw , Ktm , Kd , vars ) 
%noise , Eq. 57 

int=@(x) int noise (x ,dw, 10 ,Ktm,Kd, vars ) ; 
val=quadgk ( int , ws— gamma, ws4gamma) ; 
out=val+( vars ( 1 ) * ws*2*gamma) / (4*pi*vars(6)); 

%num 

function [ out]=num(w,dw, vars ) 

%effective occupation number with coupling , in Eq.56 

out = (2*vars (5)./ ( (w— (dw+vars (2)) — vars(3))."2+vars(5)"2)).*... 

( 2 . * n (w— (dw+vars (2)) ,vars) + l) + (2*vars (5)./ ((( dw+vars (2) ) — w— . . . 

vars(3))."2+vars (5)"2)).*(2.*n ( (dw+vars (2))— w, vars) + l); 

%params 

function [ out ] ^pararns ( 10 , dw , Ktm, Kd, vars ) 

%Numerical solution for width of signal and noise distribution 
%properties , i . e central frequency , width , . . . 
out=zeros ( length ( 10 ) , 4 ) ; 

opts=optimset ( ' MaxFunEvals ' ,1E5, 'TolX ' ,lE-25); 
for li = l:length(IO) 

[xmax, b] = fminbnd ( (@(x)— G(x , 10 (k) ,dw,Ktm,Kd, vars ) . * . . . 

num(x ,dw, vars ) ) , dw+vars (2) + vars(3)— 9000*vars (8) , dw+vars ( 2 ) + . . . 
vars(3) + 12000*vars (8) , opts ) ; 
ymax=— b ; 

[nmax , nb] = fminbnd ((@(x)— intnoise (x, dw, 10 (k ) ,Ktni, Kd, vars ) ) ,dw+ . . . 

vars (2) + vars (3) -9000* vars (8) , dw+vars (2) + ... 

vars (3) + 12000* vars (8) , opts ) ; 
nsemax=— nb ; 

limax=fzero ( (@(x)G(x , 10 (k) ,dw,Ktm,Kd, vars ).*... 

num(x ,dw, vars)— yniax/2) ,xmax, opts ) ; 
nsniax=fzero ((@(x) intnoise (x ,dw, 10 (k) ,Ktm,Kd, vars ) — . . . 

nsemax / 2) ,nmax, opts ) ; 
g=abs (xmax— hmax ) ; 
li=abs (nmax— nsmax ) ; 
out (k,l)— xmax; 
out(k,2) = g; 
out (k,3) = nmax; 
out (k,4) = h; 

end 

%s ignal 

function [ out] = signal (ws , gamma, 10 ,dw,Ktm,Kd, vars ) 
%signal , Eq.56 

int=(@(x)G(x , 10 ,dw,Ktm,Kd, vars ) . *num(x ,dw, vars ) ) ; 
out=quadgk ( int , ws— gamma, ws+gamma) ; 

%up b n d 
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function [ out]=:upbnd (dw) 

%upper boundary for bistable region, Eq.68 

out = (l/2).*dw. ' (3/2).*(l + 3.*dw."(-2) + (l-(l./dw) . 



(2)). "(3/2)). "(1/2) 



B.2 Simulation of Full Trilinear Evolution 



function [] — run_new_nic () 

%Monte— Carlo simulation of a non— degenerate parametric amplifier 
%driven by a harmonic oscillator initially in a coherent state. 
%Coupling is assumed to be to a zero temperature bath. 
%Requires Quantum Optics Toolbox 
%Paul D. Nation, Dartmouth College 

%Last modified April 2010 (code cleanup for thesis) 
N0=11; %number of basis states for zeroth mode 
Nl = ll; %number of basis states for first mode 
N2=N1 ; %number of basis states for second mode 
K=l; %value of parametric coupling strength (rate) 
gammaO = 0.0; 
gammal = 0.0; 
gainma2 = 0.0; 

alpha=sqrt ( 3 ) ; %co/ierenf state amplitude 
epsilon =0.5 i ; %squeezing parameter 
tfinal =4.0; 
dt=0.02; 

1 1 i s t =0: dt : t f i n al ; %evaluation times for evaluating differential equation 
t a u 1 i s t =K . * 1 1 i s t ; 

ntraj=l; %number of trajectories to run 
%de fining lowering operators 

aO=tensor (destroy (NO) , ident it y (Nl ), identity (N2 ) ) 
al=tensor ( identity (NO) , destroy (Nl ), identity (N2 ) ) 
a2=tensor ( identity (NO) , identity (Nl) , destroy (N2) ) 
%define quadratures 
XO=(aO+aO ')/2 
YO=(aO-aO ')/2 
Xl=(al+al ')/2 
Yl=(al-al ')/2 
X2=(a2+a2 ')/2 
Y2=(a2-a2 ')/2 i ; 
p=(al— li *a2 ) ; 

%define number oeprators for modes 0—>2 
numO=aO ' * aO 
numl=al ' * al 
num2=a2 ' * a2 
nump=p ' * p ; 

%dis sip ativ e operators for zero— temp . bath 
CO=sqrt (2*gamma0) * aO : 
Cl=sqrt (2*gammal ) * al 
C2=sqrt (2*gamma2) * a2 



109 



Appendix B: Simulation source codes 



%inital state for system: coherent state for mode and vacuum for 1&2 

vacuum^tensor ( ( basis (NO ,1)) , basis (Nl,l) , basis (N2 , 1 ) ) ; 

D=expm( alpha* aO ' — conj ( alpha )* aO ) ; %mode displacement operator 

%S=expm (0 . 5* conj ( epsilon)*a0"2~0.5*epsilon*(a0')''2); 

%D=1; 

S = l; 

init_state =D* S * vacuum ; 

%inter action picture Hamiltonian 

H=li*K*(aO*al ' * a2 '-aO ' * al*a2 ) ; 

%effective non— unitary Hamiltonian (includes losses) 
Heff=H-0.5*l i *((C0'*C0) + (C1'*C1) + (C2'*C2)) ; 

% %options for solver (needed to prevent error accumulation) 

options .lmm='BDF' ; 

options . iter='NEWroN' ; 

options . mxstep = 100000; 

options . reltol=le — 6; 

options . abstol=le — 6; 

%call to monte—carlo solver 

mc2file('test.dat', — li*Heff,{C0,Cl,C2},{},init_state ,tlist ,ntraj , options) 
mcsolve( 'test .dat ' , 'out.dat ' , 'clix .dat ' ); 
fid=fopen ( ' out . dat ' , ' rb ' ) ; %open data file 

%init . arrays for expectation values of num. oper . at tlist times 
photon_numO=zeros ( 1 , length (tlist ) ) ; 
photon_numl=zeros ( 1 , length (tlist ) ) ; 
photon_num2=zeros ( 1 , length (tlist ) ) ; 
photon_nump=zeros ( 1 , length (tlist ) ) ; 

%init matricies for prob . of nth num. state at times tlist 

PnO=zeros (NO , length (tlist ) ) ; 

Pnl=zeros (Nl, length( tlist )); 

Pn2=zeros (N2, length ( tlist )); 

%init von—Nuemann entropy 

SO=zeros ( length (tlist ) ) ; 

Sl=zeros ( length (tlist ) ) ; 

S12=zeros (length ( tlist ) ) ; 

S2=zeros ( length (tlist ) ) ; 

Sth=zeros ( length (tlist ) ) ; 

SthO=zeros ( length (tlist ) ) ; 

Sth2=zeros(length( tlist )) ; 

%initilize matricies 

Fidl=zeros ( length (tlist ) ) ; 

FidO=zeros ( length (tlist ) ) ; 

Tdist=zeros ( length (tlist ) ) ; 

TdistO=zeros(length( tlist )); 

thernial=zeros (Nl ) ; 

thernialO=zeros (NO ) ; 

variance=zeros ( 1 , length (tlist ) ) ; 

statepunip=zeros ( 1 , length (tlist ) ) ; 

s t at esi gnal=zeros ( 1 , length (tlist ) ) ; 

expecG=zeros ( 1 ,N0 ) ; 

VOX=zeros ( length (tlist ) ) ; 
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VOY=zeros ( length ( 1 1 i s t ) ) ; 
for k = l:ntraj 

if gettraj ( fid , error (' Unexpected ^data „ in „ fi 1 e ') ; end 
psi=qoread (fid , dims (init_state) ,size(tlist )); %readout state vectors 
%exp ectation values of number operator in modes 0—>2 
photon_numO=photon_numO+expect (numO , psi ) . / norm( psi ) . " 2 ; 
photon_numl=photon_numl+expect (numl , psi ) . / norm( psi ) . " 2 ; 
photon_num2=photon_num2+expect (nuni2 , psi ) . / norm( psi ) . " 2 ; 
photon_nump=photon_nunip+expect (nump , psi ) . / norm( psi ) . " 2 ; 
%calculate avg . probability of being in n—th number state at tlist 
for j =1: length ( 1 1 i s t ) 

%p artial trace over modes 

pO=p trace (psi{j},l)/trace(ptrace(psi{j } ,1)) 
pl=p trace (psi{j } ,2)/trace(ptrace(psi{j } ,2)) 
p2=p trace (psi{j } ,3)/trace(ptrace(psi{j } ,3)) 
%prob ability distributions over number state 



of each mode 



elemsO=diag ( full (pO ( 
elemsl=diag ( full (pi ( 
elems2=diag ( full (p2 ( 



))) 
))) 
))) 



%distribution 
%distribution 
%distribution 



of pump states 
of mode 1 states 
of mode 2 states 



PnO ( : , j )=PnO ( : , j ) + elemsO ( : ; 
Pnl (: , j)=Pnl (: , j) + elemsl ( : \ 
Pn2 ( : , j )=Pn2 ( : , j ) + elems2 ( : \ 
%expectation values 
nl=expect (numl , psi { j } ) ; 
nO=expect (numO, psi{ j } ) ; 
n2=expect (num2 , psi { j } ) ; 

variance (j) = ( expect (numl "2,psi{j}) — nl"2)/nl; 
[AO,BO] = eig(full(pO(: ,:))); 
bO=diag(BO) ; 
ex=zeros ( 1 ,N0 ) ; 
for q = l:NO 

for p = l:NO 

ex(q) = ex(q)+AO(p,q)*conj ( AO(p , q ) ) * ( p- 1) ; 

end 

expecO ( q)=ex ( q ) ; 

end 

%calculates log of matrix (in proper way) 
[Al,Bl] = eig(full(pl(: ,:))); 
lnO=AO*diag(log(bO))*AO'-l; 
lnl=Al*diag(log(diag(Bl)))*Al'-l; 
SO(j) = SO(j)-trace(full ( pO (:,:))* InO ) ; 
Sl(j) = Sl(j)--trace(full (pi (: ,:))* Inl ); 
S2(j)-S2(j)+(l-trace(p0*p0)); 
for i=l:Nl 

%p opulates diagonal elements of thermal matrix with <nl> 
thermal (i , i ) = (l + nl ) ~ - 1* (nl /(1 + nl ) ) ' ( i -1); 



end 
for 



i=l:NO 

%p opulates diagonal elements of thermal matrix 
thermalO (i , i ) = (l+nO) ' - 1* (nO/(l + nO ) ) ^ ( i -1); 



ith <nl> 
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end 

rho=pl ; 
rhoO=pO; 

%sqrt of density matrix 
sqrt_rhoO=(full (rhoO (: ,:)))" (1/2); 
sqrt_r ho =( f ull (rho (:,:)))" (1/2) ; 
%fidelity calculation 

Fidl (j) = Fidl (j) + trace(( sqrt_rho*thermal*sqrt_rho ) " (1 / 2)); 
FidO( j) = FidO ( j) + trace (( sqrt_rlioO*thermalO*sqrt_rhoO ) " (1/2)) ; 
xx=log(l/nl + l); 
xy=log(l/nO + l); 

%thermal entropies 

Sth(j) = Sth( j)-(log(l-exp(-xx)) + xx*(l-exp(xx))'( -1)); 
SthO(j) = Sth2(j)-(log(l-exp(-xy)) + xy*(l-exp(xy))"(-l)); 
%quadrature variance 

VOX(j)=VOX(j ) + ( expect (XO ^ 2 , psi { j })- expect (XO , psi { j } ) ' 2 ) ; 
V0Y(j)=V0Y(j) + (expect(Y0^2,psi{j})-expect(Y0,psi{j })'2); 

end 

end 

fclose(fid); %close data file (prevents errors) 
%normalize by number of traj . taken 
%avg . num. of photons in each mode at times tlist 
avg_photonO=photon_numO/ nt raj ; 
avg_photonl=photon_numl/ ntraj ; 
avg_photon2=photon_nuin2/ ntraj ; 

%prob ability of mode 0,1,2 being in the ith~^ state for the jth elem . in 
%tlist 

PnO=PnO./ ntraj ; 
Pnl=Pnl ./ntraj ; 
Pn2=Pn2./ ntraj ; 
Fidl=Fidl ./ntraj ; 
FidO=FidO ./ ntraj ; 
SO=SO./ ntraj ; 
S1=S1 ./ ntraj ; 
S12=S12. /ntraj ; 
S2=S2 ./ ntraj ; 
Sth=Sth . / ntraj ; 
Sth2=Sth2 ./ ntraj ; 
VOX=VOX./ ntraj ; 
VOY=VOY./ ntraj ; 
%quadrature claculationws 
ql=4*V0X-l; 
q2=4*V0Y-l; 

end %ends program 
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